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Abstract 

In  this  final  report  a  summary  of  our  activities  with  regard  to  HF-UHF  propagation  pre¬ 
diction  over  rough  terrain  is  provided.  Our  research  activities  in  the  HF-UHF  Propagation 
Prediction  can  be  categorized  into  three  categories:  (1)  development  of  a  theoretical  method 
to  predict  the  field  propagation  from  a  small  dipole  over  a  half-space  dielectric,  (2)  Numerical 
Wavelet-based  Method  of  Moments  approach  to  predict  the  scattering  from  random  rough 
surfaces  ,  and  (3)  Numerical  Iterative  Physical  Optics  technique  to  predict  the  scattering 
from  random  rough  surfaces  that  have  low  to  moderate  rms  slopes. 

The  accurate  prediction  of  radio  wave  coverage  at  HF-UHF  frequencies  over  irregular 
terrain  features  is  of  importance  in  the  design  and  development  of  low-cost,  low  power 
communications  system.  The  terrain  effects  consist  of  multipath,  diffraction,  scattering 
and  depolarization  of  the  electromagnetic  wave.  Current  methods  of  prediction  are  overly 
simplistic  and  tend  to  neglect  phenomena  which  have  a  significant  effect  on  radio  wave 
propagation  modeling.  We  investigated  a  number  of  techniques  in  order  to  predict  accurately 
the  propagation  over  rough  terrain  in  the  HF-UHF  range. 

1  Introduction 

With  the  rapid  expansion  of  technology  for  mobile  and  wireless  systems  an  accurate 
method  for  prediction  of  radio  wave  propagation  in  various  environments  has  become  essen¬ 
tial  in  the  design  and  development  of  efficient,  low-cost,  low-power  communications  systems 
which  can  operate  in  these  rapidly  varying  environments.  A  problem  of  particular  interest  is 
the  prediction  of  High  Frequency,  Ultra  High  Frequency  (HF,  UHF)  radio  wave  propagation 
over  irregular  terrain,  where  there  may  be  no  line  of  sight  (LOS)  path  and  the  received  signal 
may  consist  of  components  of  multipath,  diffraction,  scattering,  and  depolarization  effects 
from  various,  and  usually  numerous  terrain  features.  Realistically  all  terrain  effects  cannot 
be  accounted  for,  either  statistically  or  deterministically,  with  the  current  level  of  computer 
technology,  however  assumptions  can  be  made  which  produce  a  computational  model  of 
realistic  size  while  maintaining  an  acceptable  degree  of  accuracy. 

Several  methods  are  currently  used  in  propagation  prediction  over  irregular  terrain.  These 
consist  of  heuristic  models,  based  on  measurements,  and  electromagnetic  models.  The  irreg¬ 
ular  terrain  models  can  be  broken  down  into  two  types,  those  that  predict  scattering  from 
small-scale  roughness,  i.e.,  surfaces  whose  irregularities  are  small  compared  to  electrical 
wavelength,  and  those  that  predict  scattering  and  diffraction  from  large  scale  irregularities 
such  as  mountains  or  hills.  The  heuristic  models  are  limited  to  very  specific  physical  condi¬ 
tions  at  the  time  the  measurement  was  made  and  the  measurement  system  attributes  such 
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as  frequency,  bandwidth,  and  polarization.  Moreover  the  cost  and  time  to  perform  these 
measurements  is  prohibitive.  Current  electromagnetic  models  assume  the  propagation  be¬ 
tween  transmitter  and  receiver  is  limited  to  a  very'  narrow  Fresnel  zone  and  the  problem  can 
be  reduced  to  a  2-D  problem,  in  a  plane  between  transmitter  and  receiver,  of  propagation 
over  intervening  obstacles.  Effects  of  surface  roughness  are  handled  in  a  strictly’  empirical 
manner.  One  of  the  more  widely  used  models,  applied  to  large  scale  irregularities,  assumes 
all  obstacles  consist  of  knife  edges,  essentially  a  screen  between  transmitter  and  receiver  and 
Physical  Optics  (PO)  diffraction  methods  are  applied.  [4]  This  method  is  referred  to  in  the 
literature  as  the  knife  edge  diffraction  model.  While  this  method  is  simple  to  implement 
and  computationally  fast,  PO  methods  are  accurate  only  at  high  frequencies,  in  the  very  far 
zone  of  the  obstacle,  and  are  invalid  in  or  near  the  shadow  boundary  and  near  the  surface  of 
the  obstacle  which  obviously  does  not  reflect  many  realistic  scenarios.  Other  high-frequency 
techniques  have  also  been  applied  such  as  Geometrical  Theory  of  Diffraction  (GTD)  for 
wedge  diffraction,  but  these  techniques  tend  to  be  cumbersome  for  many  obstacles. 

Due  to  the  limitations  of  both  heuristic  and  PO  models  more  rigorous  electromagnetic 
models  have  been  investigated  in  recent  year.  Several  methods  have  been  proposed  for 
modeling  of  large  terrain  features,  including  parabolic  equation  techniques  which  assume 
cylindrical  wave  scattering  [5],  and  integral  equation  (IE)  techniques  [6]  [7].  IE  techniques 
have  been  an  area  of  major  research  in  electromagnetics  in  recent  years.  An  integral  equa¬ 
tion  is  formulated  by  enforcing  boundary  conditions  on  the  surface  of  the  scattering  object. 
The  unknown  quantities  in  the  IE  are  the  surface  currents.  These  are  solved  for  either  by 
numerical  or  iterative  techniques.  This  method  accounts  for  all  electromagnetic  phenomena 
and  interactions  and  is  accurate  to  the  degree  of  accuracy  of  the  solution  technique.  When 
a  numerical  solution  is  sought,  this  technique  is  limited  by  the  electrical  size  of  the  problem. 
While  advances  in  computer  technology  have  made  the  solution  of  larger  problems  possible 
when  using  this  technique,  the  size  of  radio  wave  propagation  problems  is  still  prohibitive 
if  all  electromagnetic  interactions  are  to  be  accounted  for.  With  this  in  mind,  techniques 
are  sought  using  IE  methods  to  reduce  the  size  of  the  radio  wave  propagation  problem  to 
a  manageable  size  while  maintaining  an  acceptable  degree  of  accuracy.  Two  techniques  to 
achieve  this  for  large  scale  irregularities  are  proposed,  a  numerical  technique  using  a  MoM 
approach  with  wavelet  basis  functions,  and  an  Iterative  PO  approach.  In  addition  a  theoreti¬ 
cal  technique  is  proposed,  using  a  perturbation  method,  to  predict  scattering  and  diffraction 
from  a  surface  with  small  scale  irregularities  when  both  transmitter  and  receiver  are  near 
the  surface. 

A  widely  used  numerical  technique  to  solve  IE  equations  is  the  Method  of  Moments 
(MoM)  technique.  This  technique  expands  the  unknown  surface  currents,  the  expansion 
consisting  of  basis  functions  with  unknown  coefficients.  The  expansion  can  be  over  the  entire 
domain  of  the  scattering  object  or  discrete  sub-domains,  where  the  scatterer  is  subdivided. 
These  sub-domains  are  at  most  1/10A  in  size  to  produced  the  required  degree  of  accuracy, 
where  A  is  the  wavelength.  The  most  widely  used  is  the  sub-domain  method.  Weighting 
functions  are  applied  to  reduce  the  average  residual  error.  Direct  application  of  MoM  pro¬ 
duces  matrices  of  order  N 2  and  number  of  computations  of  order  A3,  where  N  is  the  number 
of  unknowns  in  the  matrix  formulation.  Both  the  order  of  the  matrix  size  and  number  of 
computations  are  unacceptable  for  a  problem  the  size  of  the  radio  wave  propagation  problem. 
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Current  research  is  directed  towards  reducing  the  problem  size  while  maintaining  accuracy 
using  both  acceptable  approximations  and  computational  techniques. 

Several  methods  have  been  proposed  to  achieve  the  goal  of  reducing  the  MoM  problem 
size  in  the  area  of  radio  wave  propagation  over  hilly  terrain.  As  stated  previously  most  meth¬ 
ods  assume  a  2-D  problem,  calculating  propagation  loss  along  a  2-D  path,  and  all  methods 
cited  assume  this  unless  otherwise  noted.  Hviid  et.al  [6]  proposes  a  method  which  is  valid 
for  vertical  polarization  only.  In  this  method  the  IE  is  formulated  using  the  Magnetic  Field 
Integral  Equation  (MFIE)  and  solves  for  magnetic  currents  on  the  surface  of  the  scatterer. 
It  assumes  a  smooth,  perfectly  reflecting  smooth  surface  and  to  simplify  the  problem  does 
not  account  for  backscatter.  Brennan  et.al,  [7]  extends  this  method  by  applying  the  Fast  Far 
Field  Approximation  (FAFFA),  proposed  by  Chew  et.al.  [8]  to  reduce  the  computational 
requirements  of  the  problem.  In  this  technique  local  groups  are  defined  in  the  problem  where 
the  near  field  interaction  between  sub-domains  is  significant  and  the  IE  is  fully  enforced.  It 
is  then  assumed  that  the  distance  between  groups  is  far  enough  to  assume  plane  wave  inter¬ 
action  between  them  and  the  full  IE  need  not  be  enforced  thus  reducing  the  computational 
size  of  the  problem.  Both  methods  described  have  limitations,  including  smooth  surface 
assumptions  and  no  backscatter  in  the  first,  and  in  the  second  the  need  to  define  near  and 
far  zone.  A  method  is  sought  to  more  fully  account  for  the  electromagnetic  interactions 
while  maintaining  accuracy  and  a  computationally  efficient  problem.  With  this  in  mind  two 
techniques  are  proposed  to  solve  the  IE,  a  numerical  method  using  wavelet  basis  functions 
to  reduce  the  problem  size  while  maintaining  the  accuracy  requirement  and  an  iterative  PO 
technique,  applied  to  the  magnetic  field  integral  equation,  which  also  produces  accurate  re¬ 
sults  while  significantly  reducing  run  time.  These  techniques  are  valid  when  the  fine  details 
of  the  terrain  are  not  significant  in  the  model  as  is  the  case  at  frequencies  up  to  UHF. 

The  MoM  techniques  described  previously  use  standard  basis  functions  such  as  pulse  basis 
functions  and  point  matching  weighting  which  essentially  enforces  the  MoM  formulation 
at  the  center  of  each  sub-domain.  Alternative  basis  functions  have  been  investigated  in 
recent  years  which  produce  acceptable  accuracy  while  significantly  reducing  the  problem 
size  and  runtime.  Wavelets  are  an  adaptive  basis  function  which  has  seen  wide  application 
in  recent  years  in  both  communications  and  signal  processing  and  electromagnetic  modeling 
applications.  In  recent  years  it  has  been  shown  that  the  orthogonal  properties  of  the  wavelet 
basis  functions  in  addition  to  having  vanishing  moments  produce  a  sparse  MoM  matrix 
which  can  be  solved  using  iterative  solver  techniques  such  as  Conjugate  Gradient  (CG),  thus 
creating  a  significant  reduction  in  problem  size  and  a  significant  speedup  in  the  run  time 
of  the  problem.  [9,  10,  11,  12]  This  method  produces  accurate  results  without  the  need  to 
define  near  and  far-field  domains  and  produces  a  full  bistatic  pattern,  while  maintaining  the 
accountability  for  all  electromagnetic  interactions  inherent  in  MoM  techniques.  This  method 
has  recently  been  applied  to  the  problem  of  scattering  from  rough  surfaces  very  successfully. 

In  the  iterative  PO  technique  the  MFIE  may  be  enforced  repeatedly  to  obtain  the  so¬ 
lution,  with  the  previous  solution  substituted  for  each  successive  iteration.  The  number  of 
iterations  is  dependent  on  the  accuracy  desired.  A  significant  advantage  to  the  iterative  PO 
solution  is  the  memory  requirements  are  of  order  N  as  opposed  to  order  N2  for  standard 
MoM. 
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2  SUMMARY  OF  ACCOMPLISHMENTS 

Because  no  single  method  is  comprehensive  in  its  results  (quick  computations,  exact,  ac¬ 
curate  over  all  ranges  of  physical  parameters,  etc.)  our  techniques  have  focused  on  providing 
a  number  of  desired  qualities  including  accuracy  of  results  and  computational  speed.  To  this 
end  we  have  focused  on  three  different  techniques:  (1)  Theoretical  results  of  scattering  in 
a  dielectric  medium,  (2)  semi-exact  scattering  results  from  random  rough  surfaces  using  a 
wavelet-based  technique,  and  (3)  fast  scattering  results  using  an  iterative  Physical  Optics 
approach  for  random  rough  surfaces  with  low  to  moderate  rms  slopes.  A  summary  of  the 
work  accomplished  in  each  area  is  given  next. 

2.1  Small  Dipole  Over  a  Random  Rough  Dielectric  Surface 

With  recent  progress  in  wireless  technology  and  ever  increasing  demand  for  reliable  low 
power  wireless  systems,  the  need  for  predicting  system  performance  has  become  an  extremely 
important  step  in  the  design  of  such  systems.  At  HF  through  UHF,  the  channel  characteris¬ 
tics  such  as  attenuation,  and  multipath  fading  statistics  significantly  affect  the  performance 
of  wireless  systems.  Prediction  of  the  channel  characteristics  can  be  accomplished  using 
physics  based  propagation  models.  For  this  purpose  precise  diffraction  models  are  developed 
and  incorporated  with  accurate  terrain  data  base  to  construct  a  realistic  propagation  model. 

In  this  study  the  problem  of  electromagnetic  wave  propagation,  excited  by  a  short  dipole, 
above  a  dielectric  ground  plane  with  an  arbitrary  dielectric  profile  and  an  irregular  interface 
is  studied.  This  investigation  is  a  natural  extension  of  the  classical  Sommerfeld  problem 
with  the  exception  of  the  random  surface  irregularities  at  the  interface  between  the  two 
dielectric  media.  Assuming  that  the  interface  profile  height  variations  are  small  compared  to 
the  wavelength,  the  problem  is  formulated  as  follows.  First,  the  bistatic  scattering  of  a  plane 
wave  illuminating  the  rough  surface  is  solved  using  a  perturbation  solution  of  an  integral 
equation  for  the  induced  polarization  current.  Analytical  expressions  for  the  coherent  field 
(mean-field)  and  incoherent  scattered  power  at  an  arbitrary  observation  point,  including 
points  near  the  interface,  are  obtained.  Then  the  solutions  for  the  mean-field  and  incoherent 
scattered  power  generated  by  a  small  dipole  of  arbitrary  orientation  and  position  are  derived 
by  expanding  the  field  of  the  dipole  in  terms  of  a  continuous  spectrum  of  plane  waves  and 
using  superposition.  The  effect  of  rough  interface  on  the  surface  waves  and  the  phenomenon 
of  depolarization  caused  by  the  rough  interface  are  studied. 

In  reality,  the  interface  between  forest  canopies  and  air  is  not  flat,  hence  it  is  not  clear 
whether  the  lateral  wave  can  be  excited  and  if  it  can  how  the  surface  roughness  affects  it.  In 
this  study,  the  effect  of  roughness  of  interface  between  canopy  and  air  on  the  wave  propaga¬ 
tion  in  forested  area  is  investigated.  Also,  an  expression  for  the  mean-field  of  an  infinitesimal 
dipole  of  arbitrary  orientation  is  derived  by  obtaining  a  partial  second  order  solution  of  the 
Born  approximation  and  a  sensitivity  analysis  is  carried  out  to  demonstrate  the  variations  of 
the  mean-field  to  physical  parameters  such  as  a  effective  permittivity,  location  of  the  dipole 
and  observation  points,  and  surface  roughness.  More  details  may  be  found  in  Appendix  A. 
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2.2  Numerical  Approaches 

2.2.1  Wavelet-Based  Method  of  Moments  Approach 

The  problem  of  electromagnetic  scattering  from  rough  surfaces  has  been  the  subject 
of  intensive  investigation  over  the  past  several  decades  for  its  application  in  a  number  of 
important  remote  sensing  problems.  Radar  remote  sensing  of  the  oceans,  soil  moisture,  and 
mine  detection  using  wideband  radars  are  such  examples.  For  these  problems,  where  the 
rough  surface  is  either  the  primary  target  or  the  clutter,  the  understanding  of  interaction  of 
electromagnetic  waves  with  the  rough  surface  is  essential  for  developing  inversion  or  detection 
algorithms.  An  exact  analytical  solution  for  random  rough  surfaces  does  not  exist. 

An  alternative  approach  for  evaluating  of  the  scattered  field  and  its  statistics  for  rough 
surfaces  is  Monte  Carlo  simulation.  In  order  to  use  Monte  Carlo  simulation  for  evaluating 
the  scattering  statistics  of  rough  surfaces  more  routinely,  computationally  more  efficient  scat¬ 
tering  codes  must  be  developed.  In  this  paper  the  application  of  wavelets  as  a  basis  function 
for  the  expansion  of  induced  surface  currents  is  considered.  Traditional  method  of  moments 
(MoM)  in  conjunction  with  Galerkin’s  method  would  require  matrix  fill  computation  time 
of  the  order  of  N 2  and  matrix  inversion  computation  time  of  the  order  of  N3  (using  Gaus¬ 
sian  elimination).  It  is  well  known  that  the  solution  of  linear  system  of  equations  can  be 
obtained  far  more  efficiently  using  search  routines,  such  as  the  Conjugate  Gradient  method, 
if  the  matrix  of  the  coefficients  is  a  sparse  matrix.  In  MoM,  the  application  of  conventional 
pulse  or  rooftop  basis  and  testing  functions  would  usually  produce  full  impedance  matrices. 
Although  the  diagonal  elements  are  usually  larger  than  the  rest  of  the  elements,  the  smaller 
elements  cannot  be  arbitrarily  thresholded  without  drastically  altering  the  resulting  scatter¬ 
ing  pattern.  The  success  of  wavelet  expansion  function  in  generating  sparse  matrices  have 
been  demonstrated  for  many  circuits  and  antenna  problems  [9,  10,  11,  12].  In  the  Monte 
Carlo  simulation  of  scattering  from  rough  surfaces  the  quantities  of  interest  are  the  statistical 
parameters,  such  as  the  mean  and  variance  of  the  scattered  field,  and  therefore  it  is  expected 
that  the  overall  accuracy  be  less  sensitive  to  the  threshold  level. 

A  comparative  study  is  carried  out  to  demonstrate  the  application  of  wavelets  for  im¬ 
proving  the  computation  time  and  reducing  computational  memory  required  for  evaluating 
the  statistics  of  the  scattered  field  from  rough  surfaces  using  the  method  of  moments  in 
conjunction  with  a  Monte  Carlo  simulation.  In  specific,  Haar  and  the  first-order  B-spline 
wavelet  basis  functions  are  applied  to  the  MoM  formulation  of  two-dimensional  rough  sur¬ 
faces  in  order  to  compare  the  computation  time  and  sparsity  for  wavelets  in  the  same  family 
but  of  higher  order.  Since  the  scattering  coefficient  (the  second  moment  of  the  backscatter 
field  per  unit  area)  is  a  gentle  function  of  the  surface  parameters  and  the  radar  attributes, 
it  is  demonstrated  that  a  relatively  high  thresholding  level  can  be  applied  to  the  impedance 
matrix  which  leads  to  a  sparser  impedance  matrix  and  faster  computation  time.  It  is  also 
shown  that  applying  a  high  threshold  level  the  coefficients  of  the  high  order  wavelets  would 
increase  out  of  proportion,  however  the  effect  of  these  current  components  averages  out  when 
computing  the  scattering  coefficients. 

The  resulting  sparse  impedance  matrices  are  solved  efficiently  using  fast  search  routines 
such  as  the  conjugate  gradient  method.  A  systematic  study  is  carried  out  to  investigate  the 
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effect  of  different  threshold  levels  on  the  accuracy  versus  computing  speed  criterion.  The 
computed  scattering  coefficients  are  compared  to  previous  results  computed  using  a  conven¬ 
tional  pulse  basis  function  as  well  as  the  existing  theoretical  solutions  for  rough  surfaces. 
It  is  shown  that  wavelet  basis  functions  provide  substantial  reductions  in  both  memory 
requirements  and  computation  time.  This  procedure  is  further  discussed  in  Appendix  B. 

2.2.2  Physical  Optics  Iterative  Approach 

Development  of  numerically  efficient  Monte  Carlo-based  models  for  simulations  of  elec¬ 
tromagnetic  scattering  from  random  surfaces  has  attained  significant  prominence  over  the 
past  decade.  A  major  stumbling  block  in  this  endeavor  has  been  the  large  memory  and  com¬ 
putation  time  requirement.  This  is  because  the  size  of  the  scatter  is  large  compared  to  the 
wavelength  and  the  Monte  Carlo  simulations  demand  computation  of  the  scattering  problem 
many  times.  Iterative  methods  offer  an  alternative  approach  when  exact  solutions  are  not 
available  and  have  been  used  in  different  electromagnetic  problems.  By  construct  evaluation 
of  iterative  solutions  are  rather  straight  forward  especially  when  the  perturbation  parame¬ 
ter  is  relatively  small.  Physical  Optics  (PO)  approximation  is  known  to  provide  accurate 
approximation  for  the  induced  surface  currents  provided  that  the  local  radii  of  curvature  at 
each  point  on  the  surface  of  scatterer  is  large  and  the  surface  is  convex.  For  concave  surfaces 
and  surfaces  with  many  adjacent  humps  multiple  scattering  drastically  alter  the  standard 
PO  current.  However  these  surface  current  variations  can  be  estimated  through  an  iterative 
process.  Unlike  Method  of  Moments  (MoM)  which  requires  matrices  on  the  order  of  N 2  to 
find  the  surface  current  of  the  sample  surface  with  N  elements,  the  iterative  Physical  Optics 
method  only  requires  memory  size  of  the  order  N.  Thus,  substantial  memory  savings  are 
realized.  Also,  since  no  solver  routine  is  necessary  in  order  to  solve  for  the  surface  currents, 
as  in  the  MoM,  substantial  time  savings  are  realized  as  well. 

Iterative  methods,  such  as  iterative  PO  offer  an  alternative  approach  when  exact  solutions 
are  not  available.  Evaluation  of  iterative  solutions  are  rather  straight  forward  especially 
when  the  perturbation  parameter  is  relatively  small.  PO  approximation  is  known  to  provide 
accurate  approximation  for  the  induced  surface  currents  provided  that  the  local  radii  of 
curvature  at  each  point  on  the  surface  of  the  scatter  is  large  and  the  surface  is  convex. 
For  concave  surfaces  and  irregular  surfaces  with  many  adjacent  obstacles  multiple  scattering 
drastically  alter  the  standard  PO  current.  These  surface  current  variations  can  be  estimated 
through  an  iterative  process  producing  a  problem  size  of  order  N  as  previously  stated.  The 
computation  time  of  this  problem  then  becomes  of  order  N  which  is  significantly  less  than 
the  MoM  technique  which  is  typically  of  order  N3. 

The  application  of  iterative  Physical  Optics  (PO)  in  conjunction  with  a  Monte  Carlo 
simulation  for  characterizing  the  bistatic  scattering  coefficient  of  random  rough  surfaces 
is  examined.  The  iterative  PO  method  offers  decreased  memory  and  computation  time 
restrictions  compared  to  the  standard  numerical  methods  such  as  the  Method  of  Moments 
(MoM).  Results  from  the  iterative  PO  method  are  compared  to  the  standard  electric  field 
integral  equation  (EFIE),  the  magnetic  field  integral  equation  (MFIE)  as  well  as  the  existing 
theoretical  solutions  for  rough  surfaces.  It  is  demonstrated  that  memory  requirements  and 
computation  time  is  significantly  decreased,  even  compared  to  traditional  MoM  techniques 
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while  providing  fairly  accurate  results  for  surfaces  with  moderate  to  low  rms  slope. 

In  addition,  this  technique  is  extended  to  the  3-D  scattering  problem,  which  is  pro¬ 
hibitively  large  for  MoM  techniques  on  all  but  the  largest  and  most  powerful  computers. 
Good  agreement  is  found  to  occur  between  analytical  methods  and  the  3-D  iterative  Physi¬ 
cal  Optics  approach.  Also,  the  iterative  Physical  Optics  approach  has  been  further  simplified 
by  using  a  large  argument  approximation  for  the  kernel.  For  these  results,  even  more  time 
reduction  is  observed. 

We  investigated  the  use  of  the  iterative  Physical  Optics  method  upon  a  variety  of  surfaces 
in  order  to  find  the  approximate  region  of  validity  for  such  a  method.  The  results  obtained 
using  the  iterative  PO  method  are  also  compared  to  the  results  found  using  the  electric  field 
integral  equation  (EFIE)  with  tapered  resistive  sheets  at  the  ends  of  the  surface  samples  as 
well  as  the  magnetic  field  integral  equation  (MFIE)  for  a  horizontally  polarized  wave.  This 
procedure  is  detailed  in  Appendix  C. 
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APPENDIX  A 

Small  Dipole  over  a  Random  Rough  Dielectric  Surface 
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Abstract  -  The  accurate  prediction  of  radio  wave  coverage  at  HF-UHF  frequencies  over 
irregular  terrain  features  is  of  importance  in  the  design  and  development  of  low-cost 
and  low  power  communications  systems.  In  this  paper  the  problem  of  electromagnetic 
wave  propagation,  excited  by  a  short  dipole,  above  a  dielectric  ground  plane  with  an 
arbitrary  dielectric  profile  and  an  irregular  interface  is  studied.  This  investigation  is  a 
natural  extension  of  the  classical  Sommerfeld  problem  with  the  exception  of  the  random 
surface  irregularities  at  the  interface  between  the  two  dielectric  media.  Assuming  that 
the  interface  profile  height  variations  are  small  compared  to  the  wavelength,  the  problem 
is  formulated  as  follows.  First,  the  bistatic  scattering  of  a  plane  wave  illuminating  the 
rough  surface  is  solved  using  a  perturbation  solution  of  an  integral  equation  for  the 
induced  polarization  current.  Analytical  expressions  for  the  coherent  field  (mean-field) 
and  incoherent  scattered  power  at  an  arbitrary  observation  point,  including  points  near 
the  interface,  are  obtained.  Then  the  solutions  for  the  mean-field  and  incoherent  scattered 
power  generated  by  a  small  dipole  of  arbitrary  orientation  and  position  are  derived  by 
expanding  the  field  of  the  dipole  in  terms  of  a  continuous  spectrum  of  plane  waves 
and  using  superposition.  The  effect  of  rough  interface  on  the  surface.  waves_and.. the 
phenomenon  of  depolarization  caused  by  the  rough  interface  are  studied. 


1  Introduction 


With  the  rapid  expansion  of  technology  for  mobile  and  wireless  systems,  an  accu¬ 
rate  method  for  prediction  of  radio  wave  propagation  has  become  essential  in  the  design 
and  development  of  efficient,  low-cost,  low-power  communications  systems  .  In  many 
communication  scenarios  where  both  the  transmitter  and  receiver  are  near  the  ground, 
shadowing  and  multipath  significantly  affect  the  signal  strength  and  the  coherent  band¬ 
width  at  the  receiver.  This  is  specifically  the  case  for  the  propagation  over  irregular 
terrain.  Terrain  irregularity  so  far  as  propagation  is  concerned  can  be  categorized  into 
two  groups:  1)  large-scale  roughness,  and  2)  small-scale  roughness.  Large-scale  terrain 
irregularities  are  generally  referred  to  terrain  irregularities  large  compared  to  the  wave¬ 
length  such  as  mountains  and  hills.  Small-scale  terrain  irregularities,  on  the  other  hand, 
refers  to  surface  roughnesses  where  the  rms  height  and  slope  are  small  compared  to 
the  wavelength  (at  HF-UHF).  These  affect  the  wave  propagation  differently;  for  exam¬ 
ple,  while  large-scale  terrain  irregularities  are  the  sources  of  shadowing  and  multipath, 
small-scale  irregularities  reduce  the  ground  reflectivity  and  produces  an  incoherent  field 
component  due  to  surfaces  scattering.  Small-scale  irregularities  also  affect  the  surface 
waves  which  are  essential  when  both  the  transmit  and  receive  antennas  are  close  to  the 
air-ground  interface. 

Determination  of  field  of  a  small  dipole  over  a  half-space  dielectric  is  a  classical  prob¬ 
lem  with  a  well-known  solution  [1].  It  is  shown  that  when  both  the  transmitter  and 
receiver  are  near  the  surface,  the  contribution  from  surface  waves  is  dominant.  In  prac¬ 
tice,  the  transmit  and  receive  antennas  of  mobile  units  are  usually  very  close  (relative  to 
the  wavelength)  to  the  ground.  Existence  of  surface  roughness  may  alter  the  contribution 
of  surface  waves  drastically.  In  this  case  the  azimuthal  symmetry  of  the  problem  may 
no  longer  be  exploited,  and  the  Sommerfeld  solution  must  be  modified  significantly.  The 
surface  roughness  generates  incoherent  scattered  field  which  is  the  source  of  depolariza¬ 
tion. 

In  this  paper,  the  effect  of  slightly  rough  surfaces  on  the  radiation  of  a  short  dipole 
is  studied.  In  what  follows,  first,  a  solution  for  the  scattered  field  (including  the  near 
field)  from  a  slightly  rough  surface  illuminated  by  plane  waves  is  formulated  in  Sections 
2  and  3.  To  investigate  the  effect  of  small-scale  surface  roughness  on  surface  waves, 
ground  reflectivity,  and  the  significance  of  the  incoherent  scatter  fields,  an  analytical 
solution  based  on  a  perturbation  theory  is  proposed.  In  this  formulation  the  perturbation 
theory  is  applied  to  a  volumetric  integral  equation  for  the  induced  polarization  current 
in  the  top  rough  layer  of  the  dielectric  interface.  The  perturbation  parameter  is  the 
normalized  rms  height  of  the  rough  surface  and  an  iterative  solution  starting  from  the 
unperturbed  problem  (dielectric  half-space  with  smooth  interface)  is  obtained.  Basically, 
the  formulation  is  similar  to  what  has  recently  been  applied  to  evaluate  far-field  scattering 
from  rough  surfaces  with  inhomogeneous  profile  when  illuminated  by  a  plane  wave  [2]. 
In  Section  4  the  solution  for  dipole  excitation  is  obtained  by  expanding  the  radiated 


(icM  of  l.lir  dipole  in  l.eims  ol  <i  ronl.intioiis  sperlrmn  of  plane  waves  and  adding  llie 
solution  for  earli  plane  wave  r  olierenlly.  St.at  ist.iral  analysis  is  carried  out  analyl  irally  for 
characterizing  the  coherent  (mean)  and  incoherent  f fluctuation j  fields.  The  results  are 
comparer!  with  the  Sommerfeld  solution  and  the  depolarization  effects  are  investigated. 


2  Polarization  Current  in  a  Slightly  Rough  Surface 

As  mentioned  earlier  the  first  step  towards  evaluating  the  field  generated  by  an  ar¬ 
bitrary  dipole  above  a  ground  plane  with  rough  interface  is  to  consider  the  plane  wave 
illumination.  To  obtain  the  scattered  field,  a  perturbation  solution  to  a  volumetric  inte¬ 
gral  equation  for  the  induced  polarization  current  over  the  top  rough  layer  of  the  surface 
is  derived  using  a  procedure  similar  to  what  is  presented  in  [2],  Figure  1  shows  the  geom¬ 
etry  of  the  scattering  problem  where  a  dielectric  half-space  with  an  arbitrary  dielectric 
profile  and  rough  interface  is  illuminated  by  a  plane  wave  from  the  upper  medium.  Sup¬ 
pose  the  surface  height  variation  is  small  compared  to  the  wavelength  (A)  of  the  incident 
wave.  The  incident  wave  with  an  arbitrary  polarization  Q  can  be  written  as 

E*(f)  =  Q  e,fco*'r  , 

where  &0  =  y-  is  the  free  space  propagation  constant,  and  kl  is  the  unit  vector  along  the 
direction  of  propagation,  given  by 

k'  =  sin  6{  cos  fax  -f  sin  sin  fay  —  cos  6{Z  —  k'x  —  zk'z  . 

To  make  the  solution  tractable,  the  permittivity  of  the  top  layer  down  to  a  depth  of  d 
is  considered  to  be  uniform,  where  -d  <  mm  {surface  profile}.  Denote  the  surface  height 
profile  by  function  z  =  A/(x,  y),  where  f(x,y )  is  a  zero-mean  stationary  random  process 
with  a  known  autocorrelation  function  and  variance  1,  and  A  <<  A  is  a  small  constant 
known  as  the  perturbation  parameter.  In  the  following  derivation,  it  is  assumed  that 
the  medium  below  the  top  layer  is  stratified,  that  is,  the  relative  permittivity  is  only  a 
function  of  z. 

In  the  absence  of  the  top  homogeneous  rough  layer,  the  incident  wave  would  be 
reflected  at  the  smooth  interface  between  the  free  space  and  the  stratified  half  space  soil 
medium.  This  reflected  wave  can  be  expressed  by 

Er(f)  =  Er(0)  e,fc°*r‘r  , 

where  kT  is  the  direction  of  propagation  of  the  reflected  wave,  given  by 

kr  =  U  -  2 (z  ■  k{)z  =  kV  +  zk[  , 

and  Er(0)  denotes  the  magnitude  and  polarization  vector  of  the  reflected  wave,  which 
can  be  obtained  from 

Er(0)  =  jr„i>rf\  +  rjijii 


Q 


Hero  rv  and  r h  arc  the  Fresnel  reflection  coefficients,  and  the  horizontal  and  vertical  unit 
vectors  arc  given  by 


k3  x  z 

X  if 


v,  =  ha  x  k3  , 


(1) 


where  the  subscript  s  can  be  i  or  r  for  the  incident  and  reflected  waves,  respectively. 
In  presence  of  the  homogeneous  rough  layer,  the  incident  and  reflected  waves  induce  a 
polarization  current  within  the  top  dielectric  layer  which  is  the  source  of  the  scattered 
field.  The  polarization  current  in  terms  of  the  total  field  and  the  permittivity  of  the  layer 
is 


J(r)  =  -ik0Y0(e  -  1)E‘  , 

where  Vo  =  4-  is  the  characteristic  admittance  of  the  free  space,  and 


(2) 


E‘  =  E‘  +  Er  +  Es  . 


The  scattered  field  Es  can  in  turn  be  expressed  in  terms  of  the  polarization  current  and 
is  given  by 


Es  =  ik0Z0  J  G(r,  r')  •  J(r')  dv'  , 


(3) 


where  G(r,  r')  is  the  dyadic  Green’s  function  of  the  half-space  stratified  medium  (in  the 
absence  of  the  top  rough  layer),  and  is  given  by  [3] 

.S(r-r') 


G(r,r')  =  —  zz- 


kl 


87T2 


I  [i rhh(kz)eikzZ  +  h(-k2)e~ikzZ  h(-kz) 

i  f  ,  e«kx  (p-p')  +[rvv(k2)e,kzZ +  v(—kz)e~,kzZ]v(—k2)\e,kzZ\  if  z  <  z' 

—  /  <rkx - - - w-  r  -  a  t  * 

T  J  k*  |/i(fc2)  ]rhh{-kz)eikzZ'  +  h(k2)e~ikzZ' I 

+v(kz)  [rvv(—kz)e,kzZ>  +  v(kz)e~tkzZ']  if  z  >  z' 


(4) 

In  (4),  kz  =  y'A:2  —  A:2  —  fc2,  kj.  =  kxx  +  kyy ,  and  h(±kz)  and  v(±k2)  can  be  obtained 
from  (1)  with  ks  =  (kxx  +  kyy  ±  k2z)/k0. 

Substituting  (3)  into  (2),  the  following  integral  equation  for  the  polarization  current 
can  be  obtained: 


oo  d+Af(x\yf) 


^yJ(r)  =  -ik0Y0(Ei  +  Er)  +  k20  JJ  J  G(r,  r')  •  J(r')  dv'  . 


(5) 
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An  approximate  solution  for  llio  integral  equation  can  be  obtained  using  a  pertur¬ 
bation  technique.  The  total  polarization  current  is  expanded  in  terms  of  a  perturbation 
series  given  by 


J(r)=]Tjn(r)An  , 


(6) 


n=0 


with  the  expectation  that  lim  Jn(r)  =  0.  The  most  inner  integral  in  (5)  is  expanded  into 

n— >00 

a  power  series  in  terms  of  A  and  then  a  recursive  set  of  equations  for  J„  are  obtained. 
These  currents  are  expressed  in  term  of  their  two-dimensional  Fourier  transforms  defined 

by 


Jn(r)  =  (2^  /  rf2kl  j"(k^)e<kiP‘  • 


(7) 


After  much  algebraic  manipulation,  analytical  solution  for  the  induced  polarization  cur¬ 
rent  to  any  desired  order  is  obtained  [2].  The  expression  for  Jo(ki,z)  is  given  by 

J0(kx,  z)  =(27r)2£(kj.  -  k'x)  \j0h(z)hi  +  Jol(z)U  +  J0z(z)z]  ,  (8) 

A  A 

where  =  z  x  A,-,  and 

Mz)  =  -  i  (‘  -  1)  Got*;.*)  [<3  •  A,]  e-W  . 

MZ)  =  ~  ‘ k^+'ki,/01-1  ~  Z)  [Q  '  fl  ' 

M*)  =  - y°(£ "  *) c-(kt- z>  IQ •  *1  «-“*■'  •  (9)- 

The  parameters  used  in  these  expressions  for  the  zeroth-order  current  are  given  by 


k\z  =  k0  y/e  —  sin2  0, ,  k\ *  =  ko  sin  0{ ,  R\ 


kx  —  kl 

Kz  Klz 


fu*  ___ 


*  *£  +  *{/  s  eki  +  k[/ 


hlu  _  (~l)n(^  -  r*)  eiklzZ  +  (Rhrh  -  1) 


= 

Cvn(kp,z)  = 


72/,  (72/,  -  r/,)  -f  (72/,r/,  -  1)  e~iklzd 
(-1)"  (r„  -  FQ  +  (/2vr„  -  1)  e~ik*'* 
Ry  (Ry  -  rv)  eiklzd  +  (Rvrv  —  1)  e~iklzd 


As  before  ?’/,  and  ?■„  denotes  the  Fresnel  reflection  coefficients  of  the  half-space  medium. 
If  the  half  space  dielectric  is  homogeneous  (rx  =  72/,  and  r„  =  /2V),  the  values  of  Ck  and 
Cl  are  one.  The  expressions  for  the  first-order  currents  are  similar  to  those  of  the  zeroth 
order,  and  are  given  in  the  Appendix. 


3  Evaluation  of  Scattered  Fields 


Substituting  the  expressions  for  polarization  currents  into  (3),  and  following  a  similar 
procedure  which  resulted  in  (6)  of  [2],  the  scattered  field  to  the  Nth  order  in  A  is  obtained: 


Es(r,  k') 


w/^{/ 


G(kx;z,z')‘ 


N-l  n 


n  n  (n  +  1)!  dz'm 

n=0  m=0  v  ' 


n+1 


fln—m 


<9z,n 


(10) 


where 

6(k±;  z,  2')  =  l-{h(kz)  [r»A(-*,)ei*-*'  +  Hh)^] 

+  v(kz)  +  6(t,)e-“-']  }e“-%  if  2  >  2' 

which  is  the  Fourier  transformation  of  the  dyadic  Green’s  function. 

Equation  (10)  can  be  divided  into  two  parts: 


Es(r,  k%  k{)  =Es/(r,  kr, V)  +  Esr(r,  ks, ti)  . 


(ii) 


(12) 


_  A  A  , 

E5-*  (r,  kr,  k *)  is  the  scattered  field  to  the  zeroth  order  in  A.  Substituting  the  zeroth-order 
currents  (9)  into  (10),  the  zeroth-order  field  is  given  by 


E&(r,'k',k<)  =e 


—  Jk0kr*r 


{A(i-:)fc(-*j)40)(A;)  +  v(ki H-K)R(°HK)}  ■  Q  .  (13) 


where  Q  and  P  are  the  polarization  vectors  of  the  incident  and  scattered  fields,  respec¬ 
tively,  and 


R f\K)  =  [(i  +  RM(k;,j)  -  i]  e~2ik',d , 

R{°\K)  =  [(1  +  /&)<#(*;.  <0-1]  .  (14) 


The  zeroth-order  solution  is  equivalent  to  the  reflected  field  from  the  original  multi¬ 
layer  medium  with  a  flat  interface  ( f(x,y )  =  0),  and  the  expressions  in  (14)  give 
the  total  reflection  coefficients  at  the  air-medium  interface.  It  should  be  noted  that 
(w(— £‘), h{— klz))  =  (vi,hi)  and  (u(fc^), h(fc*))  =  (vr,hr).  The  superscript  ”/”  in  (13) 
denotes  the  flat  interface. 
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E*r(r ,  k*Jcl)  is  flic  liiglirr-ordn  scattered  field,  whic  h  only  exists  when  the  surface  is 
rough.  Substituting  the  |iolari/,;il  ion  ciirrruls  into  (ID)  mid  iifter  some  algebraic  m;mi|> 
illations,  the*  following  expressions  for  the  scattered  field  is  obtained 


-W/a 


r-r-^.r.n  =  7—  E  E 

A/  _  1  J 


t:'k (/v‘„B“,)(i/tnr 


A/  =  I 


n=0  m=0 


(/V  —  »)! 


+  (-1  r]c^(kp,d)  +  *(*,)/(*,)[*:  -  (-1  r]cvjkp,d)-i 

K0 


+i(k,)z  [n„  +  (-\r]cz+Ak„d)f 

Ka 


/  aW-n-ra-l  \  /V_n 


(15) 


Note  that  (15)  is  valid  for  all  points  z  >  d,  including  observation  points  very  close  to  the 
surface.  In  a  special  case  where  the  observation  point  is  far  from  the  surface  (2  >>  d), 
the  stationary  phase  approximation  can  be  used  for  evaluating  the  integral  in  (15).  The 
far-held  expansion  for  the  scattered  held  can  be  written  in  terms  of  a  scattering  matrix 
which  depends  on  the  direction  of  observation.  Comparing  (15)  with  (14)  and  (15)  in  [2], 
EpQ  can  be  expressed  in  terms  of  scattering  matrix  elements: 

±(A* 

N=  1  J 

sr’(k,k,)  +  b(k,)h(-K)  S<^(k,k,-)  +  «(*.)«(-*•)  S'Hk.ki)}  •  Q  , 


EZQ{r,k\k‘)  = 


2n 


A  p  •  5^>(k,k.) + h(k,)v(-ki) 


(16) 


where. the  scattering  matrix  elements  in  (16)  are  given  in  [2].  The  near-held  expression 
in  (16)  can  be  interpreted  as  the  superposition  of  the  scattered  helds  from  all  different 
directions  denoted  by  k.  The  integrand  corresponding  to  |kx|  <  &o  can  be  interpreted 
as  the  upward  propagating  waves  emanated  from  the  surface.  When  |kx|  >  ko,  the 
corresponding  waves  are  non-propagating  which  are  known  as  the  surface  wave  whose 
contributions  are  conhned  in  the  vicinity  of  the  interface.  It  should  be  emphasized  that 
the  quantities  of  interest  are  the  statistical  mean  and  the  standard  derivation  of  the 
scattered  held.  These  surface  waves  are  caused  by  the  rough-surface  scattering  and  does 
not  exist  when  the  surface  is  hat. 

Performing  the  ensemble  averaging  of  (12),  it  is  found  that,  up  to  the  second  order, 


(E*(r,  kr,ic'))  *  E^(r,  kr,k')  +  (Ear(2)(r ,fcr,fc*')>  • 


(17) 


Here  (Esr(2\v,kr,k')),  like  E5-f(r ,kr,k')  in  (13),  can  be  expressed  as 

=«“»*'•'  p ■  +  H'k:)H-K)t?(K)}  ■  Q  ,  (t«) 


t 


r 


where  R^\k'p)  and  ll\2\k')  are  also  given  in  the  Appendix.  Since  f(x,y)  is  a  zero  mean 
Gaussian  process,  it  can  he  shown  that  the  average  of  the  odd-order  fields  vanish.  The 
next  term  of  the  coherent  field  is  the  fourth-order  E4r,  which  will  be  ignored  due  to  the 
assumption  of  the  slight  roughness. 

For  the  evaluation  of  the  incoherent  scattering  power  (variance  of  the  field),  only  the 
first-order  scattered  field  is  retained.  Re-arranging  (16),  we  have 

£K,(r, *■,*•')  =  £  J (PkiWkx.kiJAFfkj.  -  ki)  .  (19) 

where  /pg(kx,k'x)  is  given  by 

Wki.ki)  =  k'AF P  '  { 5ff(k,k,)  +  *(*,)*(-*:)  sfl’fk.k.-) 

+  v(k,)h(-ki)  S<;)(k,k,)  +  t.(kI)i(-fc;)  S£>(k,kj)}  -Q  .  (20) 

Noting  that 

A2(F(kx)F-(k'±))  =  AI(F(k1)F(-ki))  =  (2*)J«(kj.  -  k'x)lV(ki)  ,  (21) 

the  incoherent  scattering  power,  up  to  the  second  order  in  A,  is  given  by 

<1^0  -  <£>«>|2>  «  (|£$f)  =  A 2  J  <?ki  |/«,(kJ.,ki±)|a  W(kL  -  ki)  .  (22) 

As  mentioned  previously,  (16)  is  expressed  as  a  continuous  spectrum  of  scattered  plane 
waves.  What  is  expressed  mathematically  by  (21)  indicated  that  these  plane  waves  are 
mutually  uncorrelated.  Therefore,  (22)  is  simply  the  integration  of  the  power  carried 
by  each  plane  wave.  For  observation  points  near  the  surface,  (22)  must  be  carried  out 
numerically  and  cannot  be  simplified  any  further.  The  convergence  of  the  integral  can  be 
examined  noting  that  W(kx)  decreases  as  |kx|  increases  and  the  fact  that  for  |kx|  >  k0, 
kz  becomes  pure  imaginary  which  causes  the  integrand  (|/pQ(kx,k^)|  )  to  decay  rapidly. 

To  demonstrate  the  effect  of  the  surface  roughness  on  the  surface  reflectivity,  a  nu¬ 
merical  example  is  considered.  Both  coherent  reflectivity  and  incoherent  reflectivity 

(<|£p«1,|  )/  I£f )  as  a  function  of  observation  point  height  are  calculated.  These  plane 
wave  illumination  examples  simulate  situation  where  the  transmitter  (receiver)  is  airborne 
and  the  receiver  (transmitter)  is  near  the  rough  interface.  Consider  a  rough  soil  surface 
with  rms  height  of  0.016m,  correlation  length  0.16m,  and  dielectric  constant  e  =  8  -f  il 
illuminated  by  a  plane  wave  generated  by  a  source  operating  at  /  =  890  MHz.  At  this 
frequency,  the  normalized  rms  height  and  correlation  length  are,  respectively,  ks  =  0.3 
and  kl  =  3.0.  Figure  2  shows  the  magnitude  of  the  zeroth  order  and  complete  second  or¬ 
der  mean-field  in  (14)  versus  incidence  angle.  In  this  simulation,  the  correlation  function 


8 


I  or  tin:  roii/>li  surface  is  assumed  Gaussian.  The  complete  second  order  solution  demon¬ 
strates  the  effect  of  surface  roughness  on  the  surface  reflectivity.  Basically,  the  surface 
roughness  reduces  the  surface  reflectivity  and  causes  a  slight  shift  in  the  Brewster  angle. 

It  should  be  noted  that  formulation  of  the  second  order  coherent  reflection  coefficients 
in  (18)  does  not  converge  for  surfaces  with  exponential  correlation  function.  This  may 
be  due  to  the  fact  that  higher  order  terms  are  excluded.  However,  this  problem  is  not 
observed  in  the  formulation  for  the  incoherent  wave.  Figure  3a  and  3b  show  compar¬ 
isons  between  the  zeroth  order  coherent  and  incoherent  reflectivities  as  the  function  of 
the  height  of  the  observation  point  for  an  exponential  correlation  function.  It  is  shown 
that  for  vertical  polarization  and  observation  point  heights  less  than  0.1  A,  the  incoherent 
reflectivity  is  significant  and  dominant  near  the  Brewster  angle.  However,  for  horizon¬ 
tal  polarization,  independent  of  the  incidence  angle,  the  incoherent  reflectivity  is  much 
smaller  than  the  coherent  reflectivity. 

Generally,  the  incoherent  reflectivities  decrease  as  incidence  angle  decreases.  This 
could  be  qualitatively  explained  by  Rayleigh  criterion  [4].  The  criterion  is  stated  as 
follows:  For  a  surface  characterized  by  a  distribution  of  irregularities  of  height  h,  if  h 
satisfies 


h  < 


A 

16  cos  6  ’ 


(23) 


where  6  is  the  incidence  angle,  the  surface  can  be  considered  smooth.  As  the  incidence 
angle  increase,  the  surface  appears  ’’more  flat”.  Therefore,  the  incoherent  scattering 
decreases.  Same  is  true  for  the  coherent  field  as  shown  in  Fig.  2  where  the  coherent 
reflectivity  approaches  unity  when  6  is  increased  to  90°. 

It  is  noticed  that  the  incoherent  reflectivities  vary  as  the  height  of  the  observation 
point  changes.  As  mentioned  previously,  the  scattered  field  can  be  decomposed  into  two 
components:  upward  propagating  waves  and  surface  waves.  When  the  height  increases, 
while  the  surface  wave  components  attenuate,  the  propagating  waves  remains  unattenu¬ 
ated.  This  phenomenon  is  demonstrated  in  Fig.  4  where  the  integrand  of  (19)  is  plotted 
in  ki  space.  The  normalized  magnitude  of  the  integrand  is  shown  in  gray  scale  over  an 
area  with  the  radius  of  2 kQ  in  the  spectral  domain.  The  propagating  waves  are  confined 
in  a  circle  of  radius  ka ,  while  the  surface  waves  are  outside  the  fc0-circle.  Figure  4(a)  and 
4(d)  show  the  incoherent  vv-  and  hh-polarized  power  spectral  densities  for  an  observation 
point  0.01  A  above  the  rough  surface  when  the  incidence  angle  is  20°.  The  integrand  is  nor¬ 
malized  with  respect  to  the  value  at  (kx,ky)  =  (ka  sin  18.8°,  0.0)  for  vv-polarization  and 
(k0  sin  21.2°,  0.0)  for  hh-polarization.  In  this  case,  most  of  the  power  is  in  the  A:0-circle, 
which  justifies  the  lack  of  sensitivity  of  the  incoherent  reflectivity  to  the  height  variation 
at  20°  shown  in  Figs.  5(a)  and  5(b).  Figures  4(b)  and  4(e)  respectively  show  the  inco¬ 
herent  vv-  and  hh-polarized  power  spectral  densities  when  the  observation  point  is  0.01  A 
above  the  rough  surface  at  incidence  angle  80°.  The  integrand  is  normalized  with  respect 
to  the  value  at  ( kx,ky )  =  (— A:o,0.0)  for  vv-polarization  and  (fco,0.0)  for  hh-polarization. 
A  significant  component  of  incoherent  scattering  is  from  the  contribution  of  the  surface 


waves  (the  area  outside  the  /^-circle.  It  is  also  noticed  that,  for  vv-polarization,  inco¬ 
herent  scattering  is  mostly  from  the  waves  around  the  backseat tering  direction.  As  the 
height  of  the  observation  point  increases  to  0.5A,  the  contribution  from  the  surface  waves 
almost  vanish,  as  shown  in  Figs.  5(c)  and  5(f). 

In  near-field  region,  all  field  components  are  present  in  general.  Decomposing  the 
field  components  into  v,  h,  and  k  components,  the  behavior  of  the  cross- polarized  in¬ 
coherent  scattered  waves  are  demonstrated  next.  Figures  5(a)-(d)  show  the  incoher¬ 
ent  cross-polarized  scattering,  including  vh,  hv,  kv,  and  kh  polarizations.  Like  the  co¬ 
polarized  scattering,  the  cross-polarized  incoherent  scattered  power  is  stronger  for  obser¬ 
vation  points  close  to  the  surface.  Also  for  v-polarized  incident  field,  the  cross-polarized 
scattering  powers  are  stronger  than  those  of  the  polarization  of  h-polarized  incident  field. 


4  Evaluation  of  Field  of  a  Short  Dipole  above  a  Rough 
Surface 


Another  problem  of  practical  importance  is  the  characterization  of  the  field  of  a 
short  dipole  above  a  rough  surface.  Consider  an  infinitesimal  current  element  given  by 
Q£(r  —  r')  where  Q  denotes  the  polarization  of  the  dipole  antenna,  and  r'  represents 
the  location  of  the  dipole.  At  the  observation  point,  r,  the  direct  radiated  field  from  the 
dipole  is  given  by  [5] 


E<i(r,  r') 


_  iZ0  f  ‘ 
47rA:ol  . 


■l  +  ik0R  +  klR2 

R 3 


Q  + 


3  -  3 ik0R  -  klR ? 


R 5 


(Q  R)R 


}■ 


(24) 


where  R  =  r  —  r'  and  R  =  |R|.  For  z  <  z' ,  (24)  can  be  expressed  in  terms  of  integral  of 
plane  waves  given  by 


,  —kZ  f  Pi^±<p-p')  r„  *  i 

E  (r,r')  =-g^T  J  - l -  [h(-kz)h(-kz)  +  v(-kz)v(-kz) j  -Qe 


-ikz(z-zf) 


(25) 


The  mean  field  can  be  obtained  by  evaluating  the  integral  of  the  mean  fields  corresponding 
to  each  individual  plane  wave  and  is  given  by 


<£Mr-r')> 


-k0Z0 


87 r2 


/ 


„-iKV 

d2k± — - — ^pg(r,  k,  I<) 


(26) 
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where  k 


kj.  -I  /•: . ,  and  K  -  kj.  -  ktz. 


Substituting  ( If)  inLo  (20),  il.  is  found  that 


(  ( r ,  r'j)  =  -!>•  KV  ,2(d  -  z’)h  H-  r')  +  f  kylk/—- 

Un*  J  k7 

o 

{o  +  H+)(%{k„d)c-2tk'dl*(kp)  +  (1  +  ltv)Cy(kp,  d)c~2ik‘d  I’^kp)  1  , 


(27) 

Equation  (27)  is  obtained  by  noting  that  the  reflection  coefficients  are  azimuthally  sym¬ 
metric  and  therefore  the  integration  with  respect  to  4>  is  carried  out  analytically  and  are 
expressed  in  terms  of  Pqikp)  and  Pq(kp)  which  are  given  in  (A. 3)  and  (A.4),  respectively. 
The  integral  in  (27)  is  a  Sommerfeld  type.  When  both  r  and  r'  are  close  to  the  ground 
and  far  from  each  other,  the  first  term  on  the  right-hand  side,  which  can  be  viewed  as 
the  negative  image  source,  dominates.  This  results  in  a  destructive  interference  with 
the  direct  wave,  and  hence  the  surface  waves,  which  is  accounted  for  in  the  integral  of 
(27),  become  dominant.  The  integral  in  (27)  can  also  be  written  in  terms  of  asymptotic 
expressions  available  in  literatures  [6].  The  numerical  technique  for  the  evaluation  of  the 
Sommerfeld  integral  is  not  discussed  here.  Interested  readers  are  referred  to  [7].  Here 
the  objective  is  to  investigate  the  significance  of  the  rough  surface  which  are  included  in 
the  integral  in  (27)  and  in  the  incoherent  scattered  field. 

The  first-order  incoherent  scattered  field  is  written  in  terms  of  superposition  of  inco¬ 
herent  scattered  fields: 


p-iK-t' 


from  which  the  incoherent  scattering  power  is  obtained  and  given  by 


(28) 


]Jw/{Km 

j  ^ki/pgfki.k'J/pgfki  +  k"  -  ki,  k"  )W(ki  -  k'J  . 


(29) 


Note  that  the  integral  in  (29)  is  six-fold  which  is  extremely  difficult  to  evaluate  numer¬ 
ically.  Practically,  the  distance  between  the  dipole  and  the  observation  point  is  large, 
which  can  be  used  to  simplify  (29).  Suppose  r  -  zz  and  |r'|  >>  A.  Using  the  stationary 
phase  approximation  to  evaluate  the  integrals  on  k'x  and  k" ,  the  incoherent  scattering 
power  can  be  obtained  from  ° 


Wkx,~ A^ofr) 
M 


W(kl+ko^:)  (30) 

In 
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Numerical  simulations  have  been  performed  to  demonstrate  the  efTect  of  the  surface 
roughness.  Consider  the  rough  surface  of  previous  example  with  parameters  c  =  8  +  il, 
ks  =  0.3,  kl  =  3.0,  and  exponential  correlation  function.  Suppose  the  infinitesimal 
current  source  is  placed  at  0.2A  above  the  surface  at  a  location  r'  =  x'x'  +  y'y'  +  0.2A z1, 
and  the  observation  point  is  at  (0,0,0.2A).  Figures  6(a)  (Q  =  P  =  z)  and  (b)  (Q  =  P  = 
y)  show  different  components  comprising  the  received  power  versus  the  radial  distance 
(%/x'2  +  y'2)  between  the  current  source  and  the  observation  point.  It  is  shown  that,  as 
the  distance  increases,  the  reflection  coefficients  approach  to  —1,  and  the  coherent  ground 
contribution  cancel  the  direct  wave.  Note  that  as  the  incidence  angle  approaches  90°,  the 
effect  of  the  surface  roughness  on  the  coherent  field  becomes  insignificant,  as  explained 
previously  when  stating  the  Rayleigh  criterion.  Under  this  circumstance,  the  dominant 
propagation  mechanism  is  the  coherent  and  incoherent  surface  waves.  H-polarized  surface 
waves  attenuate  rapidly  [8].  Therefore,  in  Fig.  6(b),  the  surface  wave  is  less  significant 
than  that  in  Fig.  6(a),  and  the  total  field  shows  obvious  interference  phenomenon  between 
the  direct  wave  and  reflected  wave.  The  incoherent  rough  surface  scattering  is  found  to 
be  somewhat  insignificant  in  both  cases,  which  was  also  predicted  in  Fig.  3  at  high 
incidence  angles. 

In  the  next  simulation,  the  distance  between  source  and  observation  point  is  fixed  at 
20A,  but  the  source  point  is  moved  on  a  circle  of  radius  20A  in  the  x-z  plane,  as  shown 
in  Fig.  7.  As  before,  the  observation  point  is  at  (0, 0,0.2A).  However  the  source  is 
at  (2OAsin0,O,O.2A  +  20A  cos  6)  with  6  €  [0,90°].  Figure  8(a)  shows  the  coherent  and 
incoherent  powers  with  Q  =  P  =  h  =  y.  When  6  approaches  90°,  the  direct  field  and 
ground  reflected  field  interfere  with  each  other  destructively,  but  the  total  coherent  field 
is  still  about  10  dB  higher  than  the  incoherent  field.  In  Figs.  8(b)  and  (c),  choosing 
Q  =  v  —  (—sin 0,0, cos 6)  whereas  P  =  x  in  Fig.  8(b)  and  P  =  z  in  Fig.  8(c).  and 
Fig.  (c)(  z)  show  the  components  of  the  received  power.  Note  that  in  both  simulations, 
the  polarizations  at  the  observation  point  are  not  suitable  to  receive  the  ground  reflected 
waves,  which  should  be  (sin0, 0,cos0).  Thus,  the  direct  field  dominates.  When  Q  and 
P  become  perpendicular  to  each  other,  the  coherent  field  diminishes,  and  the  incoherent 
field  becomes  significant. 

5  Conclusions 

The  radiation  of  a  short  electric  dipole  above  a  slightly  rough  surface  is  studied.  This 
investigation  is  a  natural  extension  of  the  classical  Sommerfeld  problem  with  the  ex¬ 
ception  of  the  random  surface  irregularities  at  the  interface  between  the  two  dielectric 
media.  In  this  paper,  the  formulation  for  the  near  scattered  field  from  a  slightly  rough 
surface  when  illuminated  by  plane  waves  is  developed  first.  A  perturbation  technique  is 
applied  to  solve  the  integral  equation  for  the  induced  polarization  current.  Analytical 
expressions  for  the  coherent  field  (mean-field)  and  incoherent  scattered  power  at  an  ar- 
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biliary  observation  point,  unhiding  points  near  the  interface,  are  obtained.  Simulations 
show  that  while  the  coherent  scattered  field  generally  dominates,  the  incoherent,  scattered 
field  is  only  essential  around  the  Mrevvster  angle  for  vv-polarizalion.  The  phenomenon  of 
depolarization  caused  by  the  incoherent  rough-surface  scattering  are  also  studied.  Gen¬ 
erally,  the  incoherent  scattered  field  become  more  significant  as  the  observation  point 
approaches  the  interface.  Then  the  solutions  for  the  mean-field  and  incoherent  scattered 
power  generated  by  a  small  dipole  of  arbitrary  orientation  and  position  are  derived  by 
expanding  the  field  of  the  dipole  in  terms  of  a  continuous  spectrum  of  plane  waves  and 
using  superposition.  Although  it  is  found  that  the  direct  and  coherent  reflected  (reflected 
+  surface  waves)  fields  are  dominant  in  most  cases,  the  incoherent  scattering  could  be 
important,  when  the  path  along  the  line  of  sight  is  obscured. 
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Appendix 


The  expression  for  the  Nth  order  polarization  currents  are  given  by 
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where  *  is  the  convolution  operator,  F(kx)  is  the  Fourier  transform  of  f(x', y'),  and  0 
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represents  n-fold  self-convolution  (0F  =  F  *  F  *•  •  •  *  F). 

The  second-order  expressions  for  reflection  coefficients  used  in  (18)  are  given  by 
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The  integrands  in  (27)  arc  given  l>y 
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Figure  1:  An  inhomogeneous  half-space  medium  with  a  rough  interface.  Left  side  of  this 
figure  shows  the  dielectric  profile. 
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Incidence  Angle  (degree) 


Figure  2:  The  magnitude  of  the  reflection  coefficients  in  (14)  versus  incidence  angle.  The 
underlying  ground  is  homogeneous,  and  the  dielectric  constant  is  e  =  8  +  tl. 
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(a)  (b) 

Figure  3:  The  comparison  between  the  coherent  and  incoherent  reflectivities  versus  in¬ 
cidence  angle.  The  underlying  ground  is  homogeneous,  and  c  =  8  +  i  1 ,  and  the  rough 
surface  has  exponential  correlation  function,  ks  =  0.3,  and  kl  =  3.0.  The  incoherent 
reflectivities  are  plotted  for  different  observation  point  heights. 
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(d)  HH,  6i  =  20°,h=0.01A  (e)  HH,  =  80°,h=0.01A  (f)  HH,  0{  =  80o,h=0.5A 


Figure  4:  The  spectral  distribution  of  the  incoherent  scattered  power  generated  from 
the  rough  surface,  (a)  and  (d)  are  respectively  for  vv  and  hh  polarization  at  incidence 
angle  20°,  and  the  height  is  0.01A.  (b)  and  (e)  respectively  for  vv  and  hh  polarization  at 
incidence  angle  80°,  and  the  height  is  0.01A.  (c)  and  (f)  are  the  same  as  (b)  and  (e),  but 
the  heigh  is  0.5A 
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(c)  (d) 

Figure  5:  The  cross-polarized  incoherent  scattering  fields  versus  incidence  angle,  (a) 
kv,  (b)  kh,  (c)  hv,  and  (d)  vh.  The  incoherent  reflectivities  are  plotted  for  different 
observation  point  heights. 
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(a)  (b) 

Figure  6:  The  electric  field  versus  the  distance  between  current  and  observation  point. 
The  underlying  ground  is  homogeneous  (e  =  8  +  il).  The  rough  surface  has  exponential 
correlation  function,  ks  =  0.3,  and  kl  =  3.0.  In  (a),  Q  =  P  =  z.  In  (b)  Q  =  P  =  y. 


Figure  7:  The  distance  between  source  and  observation  point  is  fixed  at  20A,  but  the 
source  point  is  moved  on  a  circle  of  radius  20A  in  the  x-z  plane.  The  observation  point 
is  at  (0, 0, 0.2A),  and  the  source  is  at  (20A  sin  0, 0, 0.2A  +  2OAcos0)  with  0  €  [0,90°]. 
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Figure  8:  Components  of  received  power  as  a  function  of  transmitter  height  for  a  fixed  dis¬ 
tance  (20A)  between  the  transmitter  and  receiver.  The  observation  point  is  at  (0,0, 0.2A), 
and  the  source  is  at  (20A  sin  0, 0, 0.2A  +  20A  cos  0).  In  (a)  Q  =  P  =  h  =  y,  and  in  (b)  and 
(c),  Q  =  v  =  (—  sin  0, 0,  cos  0),  and  P  =  x  and  P  =  I,  respectively. 
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Abstract 

In  this  paper,  the  problem  of  electromagnetic  wave  propagation  in  a  realis¬ 
tic  forest  environment  at  HF  through  UHF  bands  is  considered.  In  particular, 
the  effect  of  the  non-planar  interface  between  the  air  and  the  canopy,  which 
has  been  ignored  in  the  previous  models,  is  examined.  An  analytical  formula¬ 
tion  is  obtained  for  the  mean  field  when  both  the  transmitter  and  receiver  are 
within  the  foliage.  This  formulation  is  based  on  distorted  Born  approximation 
and  accounts  for  the  surface  roughness  that  exists  between  the  canopy  and  air 
interface.  It  is  shown  that  the  surface  roughness  attenuates  the  so  called  lat¬ 
eral  wave  which  is  the  dominated  source  of  the  field  at  the  receiver  locations 
far  from  the  transmitter.  It  is  also  shown  that  this  attenuation  rate  increases 
when  the  rms  height  of  the  surface  roughness  is  increased. 


1  Introduction 

With  recent  progress  in  wireless  technology  and  ever  increasing  demand  for  reliable 
low  power  wireless  systems,  the  need  for  predicting  system  performance  has  been 
become  an  extremely  important  step  in  the  design  of  such  systems.  At  HF  through 
UHF,  the  channel  characteristics  such  as  attenuation,  and  multi-path  fading  statistics 
significantly  affect  the  performance  of  wireless  systems.  Prediction  of  the  channel 
characteristics  can  be  accomplished  using  physics  based  propagation  models.  For 
this  purpose  precise  diffraction  models  must  be  developed  and  incorporated  with 
accurate  terrain  data  base  to  construct  a  realistic  propagation  model. 

Since  a  large  portion  of  earth’s  surface  is  covered  with  vegetation,  understanding 
of  electromagnetic  wave  propagation  through  foliage  is  of  great  importance  in  devel¬ 
opment  of  a  comprehensive  physics  based  propagation  model.  It  is  well  known  that 
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at  HF  through  UHF  frequencies  the  direct  propagation  between  two  distant  points 
within  a  forest  experiences  significantly  less  attenuation  than  that  predicted  by  mod¬ 
els.  That  is  the  predicted  attenuation  of  the  direct  wave  far  exceeds  the  attenuation 
values  obtained  from  experimental  results.  A  very  interesting  model  that  describes 
the  phenomena  and  attributes  the  wave  propagation  in  forest  to  lateral  wave  was  first 
developed  by  Tamir  [2].  In  this  approach,  the  forest  is  modeled  by  a  homogeneous 
half-space  dielectric  medium  with  a  planar  interface  and  a  permittivity  equal  to  the 
effective  dielectric  constant  of  the  foliage  medium.  This  effective  dielectric  constant 
can  be  obtained  using  a  dielectric  mixing  formula  [3]  or  the  formulation  for  calcula¬ 
tion  of  propagation  constant  in  random  media  [4,  5,  6].  The  former  is  appropriate  for 
low  frequencies  when  typical  dimensions  of  the  constitutive  particles  are  small  com¬ 
pared  to  the  wavelength  whereas  the  latter  is  useful  for  sparse  media  and  accounts 
for  scattering  losses.  Tamir’s  formulation  gives  the  field  of  dipole  within  the  effec¬ 
tive  medium  at  an  observation  point  near  the  interface  using  an  asymptotic  integral 
evaluation  [11,  12].  This  solution  shows  that  the  field  at  the  observation  point  is 
dominated  by  a  ray  emanated  from  the  source  point  and  traveled  in  the  direction  of 
the  critical  angle  toward  the  interface  and  along  the  interface  before  leaving  it  in  the 
direction  of  critical  angle,  in  order  to  reach  the  observation  point.  Due  to  renewed  in¬ 
terest  in  wireless  communications  at  HF  through  UHF,  recently  the  problem  of  wave 
propagation  in  forested  areas  has  gained  some  prominence.  The  origin  Tamir’s  model 
is  extended  in  [7,  8,  9]  by  representing  the  forest  by  a  two-layer  anisotropic  dielectric 
medium.  These  models  accounts  for  the  anisotropy  that  exists  in  the  trunk  layer  and 
to  a  much  lesser  extent  in  the  canopy  layer.  It  is  shown  that  upto  UHF  frequencies 
the  effect  of  anisotropy  can  be  ignored  [7]. 

In  reality,  the  interface  between  forest  canopies  and  air  is  not  flat,  hence  it  is  not 
clear  whether  the  lateral  wave  can  be  excited  and  of  it  can  how  the  surface  rough¬ 
ness  affects  it.  In  this  paper,  the  effect  of  roughness  of  interface  between  canopy  and 
air  on  the  wave  propagation  in  forested  area  is  investigated.  In  section  2,  analytical 
formulation  using  the  volumetric  integral  equation  in  conjunction  with  the  distorted 
Born  approximation  is  presented.  An  expression  for  the  mean-field  of  an  infinitesimal 
dipole  of  arbitrary  orientation  is  derived  by  obtaining  a  partial  second  order  solu¬ 
tion  of  the  Born  approximation.  In  section3,  a  sensitivity  analysis  is  carried  out  to 
demonstrate  the  variations  of  the  mean-field  to  physical  parameters  such  as  a  effective 
permittivity,  location  of  the  dipole  and  observation  points,  and  surface  roughness. 
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2  Analytical  Formulation 

In  this  section,  an  analytical  formulation  for  the  evaluation  of  mean  field  generated 
by  a  small  dipole  in  a  homogeneous  half-space  dielect  ric  medium  with  rough  interface 
is  derived.  The  effective  permittivity  of  forest  canopies  is  slightly  higher  than  that  of 
free  space  since  the  vegetation  particle  density  is  very  small  [2].  In  this  case,  distorted 
Born  approximation  may  be  used  to  evaluate  the  scattered  fields.  The  geometry  of 
the  diffraction  problem  is  shown  in  Fig.l  where  the  dipole  and  the  observation  point, 
are  respectively  located  at  h  and  h'  inside  a  canopy  with  effective  dielectric  constant 
Z\.  The  envelope  of  the  canopy-air  interface  is  denoted  by  a  random  process,  d(x,y). 
The  permittivity  of  the  upper  medium(air)  is  denoted  by  £2 ■  We  modify  the  problem 
by  extending  the  canopy  to  z  =  0  plane  and  assume  that  in  this  region,  there  exists  a 
volumetric  polarization  current  J  =  ik\Y\{ei  —  £i)E',  where  E'  =  E!  +  Er  +  Es.  The 
unknown  total  internal  field  can  be  obtained  from  the  following  integral  equation 

/oo  poo  pO  _ 

/  /  G(Tf,7)Et{r')dx,dy'dz' 

■oo  J  —  00  J —d(x,y) 

where  E*  is  the  incident  field,  Er  is  the  reflected  wave  from_the  planar  interface,  z 
=  0  and,  Es  is  the  scattered  field  generated  by  J  itself  and  G  is  the  dyadic  Green’s 
function  of  the  half  space  dielectric  medium(see(2)  [1,  10]). 

To  the  first  order  in  (£2  —  £1),  the  Born  approximation  can  be  used  to  the  scattered 
field 

/°o  poo  pO 

/  /  0{r, 7)  ■[&{?) +  Er(r'))dx'dy'dz'  (1) 

■00  J  —00  J  —d(x,y) 

In  (1),  d(x,y)  is  a  two  dimensional  random  process  describing  the  interface  between 
the  canopy  and  air  and  is  assumed  to  be  Gaussian  with  a  mean  value  of  m(a  positive 
number)  and  a  standard  deviation  of  a.  Noting  that  the  z-dependence  of  E!  and 

Er  are  of  the  form  of  exp(±jk2zz'),  where  k2z  =  \jk\  —  klx2  —  kly2,  the  integration 
with  respect  to  z'  can  be  carried  out  analytically.  This  integration  would  result  in 

an  algebraic  expression  in  terms  of  d(x,  y)  which  is  amenable  for  the  calculation  of 

— * 

statistical  average  of  Es.  It  is  shown  that  distorted  Born  approximation  provides 
a  more  accurate  solution  for  Es.  In  this  approximation,  a  phase  correction  term  is 
incorporated  into  the  approximate  expressions  for  the  internal  field  to  account  for 
the  difference  in  the  propagation  constant  between  the  air  and  the  canopy.  The 
explicit  expressions  for  the  incident  and  reflected  fields,  including  the  appropriate 
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(2) 

(3) 


phase  correction  terms  for  a  plane  wave  illumination,  are  given  by 

E'  =  P(k[  »'+*!*']  ■  e,'(fcL-fcU)m 

Er  =  P{-k\.)  RvM  ei[k*x'+k'y  y'-k*~]  ■  ei{k^-k^)m 

where  P(/liz)  is  a  unit  vector  indicating  polarization  of  the  incident  wave  (P  =  v  or 
h)  and  Rv>h  is  the  Fresnel  reflection  coefficient  for  vertical  or  horizontal  polarization. 
The  unit  horizontal  and  vertical  polarization  vectors  are  defined  with  respect  to  the 
plane  of  incident  and  given  by 

7  t  A 

h{klz)  =  -A  —  Z  ,  h{-ku)  =  h(kiz) 

|  k\  x  z  | 

v(klz)  =  h{klz )  x  k[,  v(-klz)  =  h(ku )  x  k[, 


where  k[  =  k\  —  2 (k\  x  z)z  denotes  the  direction  of  the  reflected  wave. 

Since  the  layer  thickness  is  not  uniform,  the  phase  correction  terms  in  (2)  and  (3) 
are  chosen  for  a  layer  with  a  uniform  thickness  m.  The  dyadic  Green’s  function  for 
the  half-space  dielectric  medium  when  2  <  ^'(observation  point  inside  the  foliage)  is 
given  by  [1,  10] 


G  =  - 


zz 


J(r-P) 


+ 


ko 2  ‘  87T2 


poo  poo  1 

—  y  ^  y  dkxdkye*k*i*-*‘ { h(-klz )  [Rhh(klzy 


e~ikl*z'  +  h(-hz)eikl*z']  +  v{-ku)  [Rvv{klz)e~ik'’z'  +  v{-klz)eiki*z' ]  }  e~ik'‘ 


(4) 


After  substituting  (4)  into  (1),  the  integration  with  respect  to  z'  is  carried  out  an¬ 
alytically.  After  performing  this  integration,  evaluation  of  (Es)  is  attempted  which 
requires  computation  of  the  term  like  (e±jsd).  Assuming  a  Gaussian  p.d.f.  for  d,  it 
can  easily  be  shown  that 


(eisd)  =  ^  e  V  erfc(-j£ |)  +  e 


(■ e~isd ) 


2e*2 

~i~erfc(-j^)  +  c“ 


'erfc(j^)\ 

*~erfc(j$ §)] 


(5) 

(6) 


where  complex  variable,  s  can  be  either  =  k2z  +  k[z  or  s2  =  k^z  —  k\z,  and  erfc(-) 
is  the  complementary  error  function  defined  by  [13]: 


erfc(z )  =  1  — 


et2+j2atdt 
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where  z  —  a  +  j3. 

Noting  that  erfc[jx )  =  1  —  ^ j  J_rei  dt  for  real  x,  it  can  readily  be  shown  that 

recJUCes  to  the  more  familiar  expression,  (e~^sd)  =  e~Jsm  3~. 

The  integration  with  respect  to  x'  and  y'  can  now  be  carried  out  analytically  which 
results  in  8(kx  —  ktx)-8(ky  —  ky).  This,  in  turn,  simplifies  the  integration  with  respect  to 
kx  and  ky.  The  final  result  is  a  plane  wave  propagating  along  k\.  Hence  the  total  field 
in  medium  2  is  the  sum  of  two  plane  waves:  1)  the  reflected  wave  at  the  hypothetical 
planar  interface  at  z  =  0,  and  2)  the  mean  scattered  field.  Hence  to  the  first  order  in 
fco2(e2-£i), 

E*  =  Er  +  (Es)  «  {Rref  +  R{1  ■>rn)ci[fci,*+*i,»-*i,«l  (7) 


where  Rrej  is  the  Fresnel  reflection  coefficient  for  the  canopy-air  boundary. 

After  some  algebraic  manipulations,  the  first  order  reflection  coefficients  (referenced 
to  z  =  0  plane)  for  the  mean  field  are  obtained  and  given  by 


d(i) 

Born (v) 


%£l)  { -(1  -  Rl  +  Rl(eiSid)  -  (e-is'd))  •  kp2  ku 

ZKlz  £i 


+—((ets*d)  -  (e-is>d))k02\ 

52  J 


,tS2m 


(8) 


p(i) 

^Bornih.) 


ko2(£2  —  £i)  f  1 


2  k[z 
Rh  , 


^j-(l-Rl  +  Rl(e^d)-(e-^d)) 


-\-—((eis2<i)  -  (e~is2d)) 
5  2 


0lS2Tn 


(9) 


Close  examination  of  Rglrn  reveals  that  Rre/  +  Rsom  is  accurate  enough  for  horizontal 
polarization,  but  it  cannot  accurately  predict  the  Brewster  angle  for  the  vertical 
polarization.  To  rectify  this  difficulty,  higher  order  solutions  must  be  obtained.  This  is 
accomplished  by  obtaining  a  partial  second  order  solution  for  the  vertical  polarization. 
The  second  order  solution  for  Es  is  given  by 

Es(2)(r)  =  k40(e2  -  £l)2  jG(r,  7)  •  jjf  &(?,?)  •  [E*'^)  +  Er(F,)]|  dv'dv  (10) 


Using  only  term  of  the  dyadic  Green’s  function  (see[4])  in  the  integrand  of 


the  inner  integral,  a  partial  second  order  solution  is  obtained. 

,(e2  —  £i)2 


Es(1'5)(r)  =  -k02 


£2 


/  G(r, 

J  V 


r)  •  H-(E‘(0+Er(F')) 


dv' 


(11) 
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The  complexity  of  this  integral  is  of  the  same  order  as  that  of  (1)  and  the  ensemble 
average  of  E4*1-5)  can  be  obtained  in  a  manner  similar  to  what  was  used  in  computation 
of  (Ed1').  Including  the  partial  second-order  term,  the  reflection  coefficient  for  the 
vertical  polarization  is  modified  as  follows. 


od-s) 

lBorn(v) 


((e“’J>  -  <e-”’d»  ■  ( 

S2  £ 2 


“ kp 2  k\z  ) 


Rl  +  Rl(e“'d)  -  (e~‘aid))  ■  ( 


-isld\ 


£ 2 


kp2 


-  *42)} e!S2m 


(12) 


2  ■  2  •  2 

where  kp  =  k'x  +  k'y  .  Now  it  can  easily  be  shown  that  at  Brewster  angle  where 
Rv  =  0,  R(Born(v)  vanishes  also. 

Having  found  an  approximate  analytical  solution  for  a  plane  wave  excitation,  the 
solution  for  an  infinitesimal  dipole  embedded  within  the  foliage  can  be  obtained  by 
expanding  the  field  of  the  dipole  in  terms  of  a  continuous  angular  spectrum  of  plane 
waves.  Without  loss  of  generality,  consider  a  dipole,  whose  orientation  is  denoted  by 
a  unit  vector  /,  located  at  ro  =  —(h  +  m)z.  The  field  of  this  dipole  within  a  uniform 
medium  with  permittivity  £j  can  be  expressed  by  [10] 


Ed(r)  = 


IkoZo  r  f°°  f  1  -  fcjjjj. 

8tt2  y.oo  J .00  [  k\z  kj 


(13) 


where  I  is  the  amplitude  of  the  sinusoidal  current  carried  by  the  dipole.  The  integrand 
of  (13)  can  be  regarded  as  E0e,ki'if  where 


IkoZo  1 

“8^~  Ml 


(k\  ■  1)9 
k\ 


is  the  amplitude  of  each  plane  wave  propagating  along  k\  =  k\/ki.  Using  superpo¬ 
sition,  the  mean  scattered  field  in  the  presence  of  the  upper  free-space  medium  with 
rough  interface,  can  be  computed  from  the  coherent  sum  of  all  reflected  plane  waves. 
That  is 


<E‘>  «  Ere/  - 


IkoZo 

87T2 


nOO  1 

TT  0 ■  +  (l  ■  )MK,)R 

■00  Kl2  L 


(1) 

Bh 


.  e-ik\z(z-m-h)ei(kix+k'yy)dkidki 


(14) 
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where  Ere/  is  the  field  of  the  dipole  in  the  presence  of  the  upper  free-space  medium 
with  a  planar  interface  at  z  =  0.  The  expression  for  Ere/  is  similar  to  the  integrand  in 
(14)  with  the  exception  that  Rb  is  replaced  with  Rrej.  Using  the  change  of  variables 


x  =  pcos(f,  y  =  ps'mp 
kx  =  kp  cos  7,  ky  =  kp  sin  7 

and  making  the  use  of  the  following  identities 

27T 

cos  mueik»pco<v-uUu  =  2mm  cos  nupjm{kpp) 

27T 

smmi/e,kpPCOS^v~^du  =  2nim  sin  rrnpjm(kpp), 
the  integral  in  (14)  can  be  simplified  to 


(Es)  = 


IkoZo  r°„  kpe-ik'^-m-h) 


47T 


/ 


dk0 


ku 


pi1) 

nBh 


{(Mkpp)  +  Ji{kpp)  cos  2 <p)lx 


l1-5)  p 

+  J2{kPp)  sm2<ply}  +  -jj-{ — Y^kpp)  ~  MkpP) cos  2 <p)lx 


k 2 

+~Y’J2(kPp)  sin  2 iply  —  ikpk\z  cos  <pJi(kpp)lz} 


A  1 

x  + 


pU) 

-ttWKp)- 


p(15)  k 2 

sin 2 (plx  +  (Jo(kpp)  -  J2(kpp) cos  2<p)ly}  +  -Jk[-{-lLj2(kpp)  • 
k2 

sin  2 <plx - 7^(Jo(kpp)  +  J2(kpp)  cos  2(p)ly  -  ikpku  sin  <pJi(kpp)lz} 

p(15) 

1  nBv 


1^ 2  [i^p^lar(cOS  (plx  “1"  sin  “h 


(15) 


Obviously  a  similar  expression  can  be  obtained  for  Ere/.  As  will  be  shown  later,  the 
accuracy  of  the  distorted  Born  approximation  degrades  as  the  mean  height  m  and/or 
the  normalized  permittivity  difference  increase.  To  rectify  this  problem  to  some 
extant,  we  use  the  fact  that  an  exact  solution  is  available  for  a  planar  dielectric 
interface.  In  fact  for  a  plane  interface^  =  0)  at  z  =  -m,  for  the  reflection  coefficient 
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of  each  plane  we  have 


Rexact(a  =  0)  =  Rrefei2k 
Hence,  for  a  rough  interface,  we  may  write 

Rexa ct{&)  ~  Rexact(,&  ~~  0)  R}3orn(&  =  0)  T  RBorn{^)  (16) 

and  substituting  (16)  into  (15),  the  modified  solution  is  given  by: 

<EL«,(<7)>  »  =  0)  -  E’B„Ja  =  0)  +  <EJ,or»)  (17) 


3  Asymptotic  Evaluation 


In  the  previous  section,  an  closed  form  solution  for  the  mean-field  of  a  dipole  inside 
a  dielectric  medium  with  rough  interface  was  obtained.  When  the  radial  distance 
between  the  observation  point  and  the  source  point  (p)  is  large,  the  integrand  becomes 
highly  oscillatory  and  therefore  accurate  numerical  evaluation  of  the  integral  becomes 
extremely  inefficient.  This  is  especially  true  when  both  points  are  near  the  interface. 
The  standard  approach  is  to  change  the  contour  of  the  integral  of  integration  by 
first  extending  the  limit  of  the  integral  over  the  entire  real  axis(— 00,00)  and  then 
using  a  change  of  variable  kp  =  k\  sin  w.  An  approximate  analytical  expansion  can 
be  obtained  by  applying  the  standard  technique  of  the  steepest  descent  [11,  12],  The 
expressions  for  the  reflection  coefficients,  in  the  tu-plane,  take  the  following  forms 


Rh 


cos  w  —  yfn  —  sin2  w  n 

/  .  ==>  My  = 

cos  w  +  v  K  —  sin  w 


k  cos  w  —  v k  —  sin2  w 
k  cos  w  +  y/K  —  sin2  w 


where  k  =  ^  <  1.  These  introduce  a  branch  cut  (with  the  branch  point  at  = 
sin-1  y/K.)  in  the  to-plane  which  may  be  encountered  by  the  steepest  descent  path. 
In  this  case,  diffracted  field  is  dominated  by  the  saddle  point  and  the  branch  cut 
contribution.  Hence,  in  general,  each  component  of  the  diffracted  field  may  be  written 
as 


~  Is  +  U(0S  —  0c)Ibc  (18) 

where  Is  and  Ibc  represent  the  saddle  point  and  the  branch  cut  contributions  respec¬ 
tively  and  U(-)  is  the  Heaviside  step  function,  9S  is  the  saddle  point,  and  0C  is  given 
by 

9C  —  Re{tU6}  —  cos-1  [sech(Im{tni,})] 
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The  saddle  point  contribution  can  be  obtained  rather  easily.  The  saddle  points  cor¬ 
responding  to  the  integrand  of  (14)  is  ivs  =  0li2  where  sin #1.2  =  costh  = 

cos  02  =  !±±h-±2in.'  ri  =  4.  (h  -f  h')2,  and  r2  =  -y/p2  +  (/;  +/?'  +  2m)2.  The  saddle 

point  contribution  is  found  to  be 

E-ad  - j  ^°Z°  [£ — . —  +  /?t,cos2^1att,+  cos<psin20i/?vT}  + 

"  e4L  (19) 

- {  RBhdih  +  COS2  02a,„+  cos  <p  sin  202#Bv/-}] 

r2 

for  subscript  i  =  x  or  y,  and 

axb  =  (1  —  cos  2^p)lx  —  sin  2(p/y,  =  —  (1  +  sin2<p)/x  —  sin2cp/j,, 

ayh  =  —  sin2<p/r  +  (1 -(- cos2<p)/,,,  =  —  sin2<p/r  —  (1  —  cos2<p)/y. 

The  z  component  of  the  field  is  given  by 


gsad  ^  [Rv  sin^i - {cos  iflx  -f  sin<pL  +  sin#i/2)  + 

4  7T  *’ 


Rbv  sin  &2 - {cos  iplx  +  sin  c ply  +  sin  02/2}] 

r2 

In  (18)  and  (19),  the  Fresnel  ( Rv,Rh )  and  Born  (Bb^Rb^)  reflection  coefficients  are 
evaluated  at  the  saddle  point.  Evaluation  of  the  branch  cut  is  far  more  complex.  In 
this  case,  the  integral  is  expanded  near  the  branch  cut  and  the  contour  is  deformed 
to  the  steepest  descent  at  the  branch  cut  using  a  change  of  variables  [11] 

cos(u>  —  0)  =  cos(wb  —  0)  +  js2 

After  much  algebra,  the  branch  cut  contribution  of  (14)  is  given  by 


?b ra  IZl  /  2j  sin  Wb 


[U(0i  -  9c)fi  {R!haih  +  R'vcos2wbaiv+  cos<pcos2wbR'Jz}  + 

U(0 2  -  0c)f2  | R'Bhaih  +  Rg'^  cos 2wbaiv+  cos  <p  cos  2wbR^  lz 

(21) 


for  the  x,  and  y  components,  and 


IZi  2  j  . 

8; r  V  P  81 


sin3'2®,  [(7(0,  -  ),)/,  R’„+  U(02  -  ec)f2R%  f 
{cos  wb  (cos  (flx  +  sin  <ply)  —  sin  wblz} 
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for  the  z  component.  In  (21)  and  (22).  ft  =  -  J<~ir' [°,lH  8,1 .  an(l  [{'  ft'  ,  and 

sin(u/b-0,  )rf  1  '  DX 

R'Bh  are  coefficients  of  the  linear  term  in  the  Taylor  expansion  of  the  each  reflection 
coefficient  at  the  branch  point  given  by 


where 


and 


A 

C 

E 

a 


K  = 


2a 


>  K  = 


2a 


k  cos  Wb  ?  COS  11/ 5 


^ Bv  2 T[k[  cos  tul  ( 9v\  +  9v 2)5  Rsh  —  2*1  cos  wl  (9hl  +  9hl ) 


9v  1 

9v  2 

ghi 

9hi 


A  Da 
k\  cos  Wb 
AEa 
ki  cos Wb 
Aa 

ki  cos  Wb 
Aa 

k\  cos  Wb 


2^C  j2k1me-j2kimcosw»  +  CB 


K  COS  Wb 


Aa 


j2k\m  — 
BC 


4  -kC  CB 


Aa 
j2k\  m 


k  cos  Wb  Aa 
-j2k1me-j2kimcosWb 
4  +  C  BC 


C 


COS  Wb 


cos  Wb  Aa 


cr2a 


R e{h{wb)  -  1},  B  —  Re{h(wb)-^-ki  cos  wb  -  j-^L}, 

2  V27T 

1  -  e-j2k'mcosWb,  D  =  k2  cos2  wb  +  S  sin2  wb, 

K 

—k\  cos 2  wb  +  —  sin2  wb,  h(wb )  =  cos2  Wbl2erfc{-j  ^  ), 

/c 

j  sin  2iCf, 
sin(n;()  —  0) 


The  above  formulation  is  valid  when  the  observation  point  6  is  away  from  the  complex 
critical  angle.  It  can  easily  be  shown  that  this  formulation  reduces  to  the  flat  boundary 
case  by  letting  <r  go  to  zero.  In  this  case,  h(wb)  =  1  which  reduces  A  and  consequently 

Rg'y^  and  R'Bh  to  zero.  The  remaining  terms  in  (21)  and  (22)  would  be  the  branch 
cut  contribution  from  a  flat  boundary.  To  extend  the  valid  region  to  observation  point 
near  the  critical  point,  (18)  may  be  modified  as 


E;  ~  I,  +  u(6  -  ec)IBC  ■  F{P) 


(23) 
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where  F(/3)  is  a  correction  factor.  This  function  can  be  obtained  by  retaining  higher 
order  terms  of  the  Taylor  series  expansion  of  the  integrand  near  the  branch  cut.  For  a 
dielectric  with  fiat  interface,  F(/3 )  is  represented  by  the  Weber  or  parabolic  cylinder 
function  as 

F(0)  =  e;(02/2+3*/8)(2£2)3/4  D_3/2(/?  +  j7?) 

where  / 3  =  2\Jklri  sin  61  ~gc .  Here  D^/2{(3  +  j/3)  is  Weber  function  of  order  -3/2  and 
is  defined  as  [11,  12], 

o^a72/4  roo 

°n(x)  = 

'  (~  n)  Jo 

As  13  increases,  F{(3)  approaches  unity  and  (23)  reduces  to  (18).  It  is,  however, 
extremely  difficult  to  modify  the  formulation  for  a  rough  boundary  in  such  simple 
manner  as  the  higher  order  derivative  of  the  special  functions  in  the  integrand  are 
difficult  to  evaluate.  The  asymptotic  expressions  clearly  indicate  that  the  diffracted 
field  due  to  the  branch  cut  contribution  (see(21)  and  (22))  decays  as  1/p2  whereas  the 
saddle  point  contribution  decays  exponentially  (ejfcir/r).  Hence  for  far  observation 
points  the  lateral  wave  is  the  dominant  source  of  the  diffracted  field.  As  will  be  shown 
later,  the  lateral  wave  for  rough  boundaries  are  weaker  (depending  on  k0a)  than  the 
lateral  wave  for  a  flat  boundary. 


4  Numerical  Simulation 

In  this  section,  numerical  examples  are  considered  to  examine  the  validity  of  the  so¬ 
lution  based  on  distorted  Born  approximation  and  to  demonstrate  the  sensitivity  of 
the  field  intensity  at  the  observation  point  to  the  canopy  parameters  such  as  effective 
dielectric  constant,  canopy-air  interface  roughness,  and  transmitter  and  receiver  posi¬ 
tions.  In  the  following  simulations,  the  transmitter  dipole  is  assumed  to  be  operating 
at  f=30MHz  having  II  =  1.  To  examine  the  effect  of  transmitter  polarization,  a  verti¬ 
cal  dipole  (/  =  z)  and  a  horizontal  dipole  (/  =  y)  are  considered  in  a  canopy  with  £1  = 
1.03+j 0.006,  and  flat  interface  at  h= 3m  (0.3A).  Figure  2(a)  and  2(b)  show  the  three 
components  of  the  diffracted  electric  field  (excluding  the  direct  field  contribution  from 
the  dipole)  as  a  function  of  distance  for  a  observation  point  3m  below  the  interface  ( h ' 
=  3m)  for  the  vertical  and  horizontal  dipole  respectively.  Also  shown  in  these  figures 
are  the  results  obtained  from  the  asymptotic  evaluation(see  (18-22)  for  k0cr  =  0). 
From  these  figures,  it  is  obvious  that  the  field  of  a  vertical  dipole  experiences  far  less 
propagation  loss  than  the  field  of  a  horizontal  dipole. 
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Figure  3  shows  the  vertical  component  of  the  total  field  (diffracted  plus  the  direct 
field)  of  a  vertical  dipole  as  a  function  of  distance  when  h  =  2m  (0.2A)  and  h'  —  lm 
(0.1A)  for  three  different  values  of  e\.  In  this  simulation,  the  imaginary  part  of  ty 
is  kept  constant  and  the  real  part  is  changed  from  1.01  to  1.05.  Figure  4  shows  a 
simulation  similar  to  that  of  Fig.3  with  the  exception  that  here  the  real  part  of  the 
effective  permittivity  is  kept  constant  and  the  imaginary  part  is  increased.  The  total 
field  calculation  for  Figs.  3  and  4  is  accomplished  using  a  numerical  integration  for 
the  observation  points  as  far  as  2km  and  for  the  further  point  the  asymptotic  method 
is  used.  As  expected,  the  propagation  loss  increases  with  increasing  the  imaginary 
part  of  e\.  Next  we  examine  the  effect  of  location  of  transmitter  and  receiver  in  the 
canopy.  It  should  be  noted  that  the  dependence  of  the  field  on  h  and  h!  is  of  the  form, 
h  +  h'.  Figure  5  show  the  dependence  of  the  field  of  a  vertical  dipole  as  a  function  of 
ko(h  +  h')  for  three  different  permittivity  values  at  fixed  lateral  distance  ( p  =  2.8km). 
Unless  the  observation  and  source  points  are  very  close  to  the  surface,  the  field  decays 
exponentially  as  a  function  of  k0(h  +  hr)  with  an  exponential  factor  proportional  to 
Im[\/1  —  «]. 

To  examine  the  accuracy  of  the  distorted  Born  solution,  a  vegetation  layer  with 
smooth  boundary  is  considered  (<x  =  0,  and  a  non-zero  m).  For  this  case  an  exact 
solution  exists.  First  the  Fresnel  reflection  coefficients  are  considered.  Figure  6(a)  and 
6(b)  compare  the  magnitude  and  phase  of  the  phase  compensated  Fresnel  reflection 
coefficient(f?/le2jfclzni)  of  a  dielectric  interface  with  ei  =  1.03  +  j 0.005  at  m  =  0.6A 
with  those  computed  by  the  approximate  distorted  Born  solution  for  perpendicular 
polarization.  The  computations  are  performed  at  30MHz  for  two  different  values  of  £\ 
and  two  different  layer  thickness  values  as  a  function  of  sin  0,-  =  kp/k0-  Similar  results 
for  parallel  polarization  are  shown  in  Fig. 7(a)  and  7(b).  Their  errors  in  magnitude 
and  phase  do  not  exceed  5%  and  5°  respectively.  The  accuracy  of  the  distorted  Born 
approximation  degrades  as  the  dielectric  constant  and  layer  thickness(m)  increases. 
To  determine  the  region  of  validity  of  this  solution,  the  exact  solution  was  compared 
with  the  distorted  Born  approximate  solution  for  a  wide  range  of  £\  and  m.  It  was 
found  that  the  accuracy  of  the  approximate  solution  improves  as  the  imaginary  part 
of  the  dielectric  constant  increases.  To  specialize  the  region  of  validity  to  the  problem 
at  hand  and  lower  the  number  of  independent  variable,  we  considered  a  canopy  with 
vegetation  particle  permittivity  ev  =  50  +  j25  and  used  Polder  Van  Santer  dielectric 
mixing  formula  [3]  to  compute  the  effective  dielectric  constant  £\  from 

£t-l  ,  _  £i  -  1 
£v  +  2fii  3(?i 

where  /  is  the  volume  fraction.  Tolerating  5%  error  in  magnitude  and  5°  error  in 
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phase  of  the  reflection  coefficient,  the  region  under  curve  shown  in  Fig.S  denotes 
values  of  m  and  £1  where  the  distorted  Born  approximation  produces  valid  results 
for  the  reflection  coefficients.  With  a  confidence  in  the  formulation  of  the  distorted 
Born  approximation,  the  effect  of  surface  roughness  on  the  reflection  coefficient  is 
considered  next.  Figure  9(a)(perpendicular  polarization)  and  9(b)(parallel  polariza¬ 
tion)  show  the  magnitude  of  reflection  coefficients  as  a  function  of  kp/k0  =  sin0,- 
and  two  different  values  of  surface  rms  height  a  =  7r/10  and  a  —  7r/5.  In  these 
simulations,  £\  =  1.03  +  j‘0.015  and  m  was  chosen  to  be  2a.  It  is  shown  that  the 
reflection  coefficient  is  drastically  reduced  by  the  surface  roughness  for  low  incident 
angle  and  for  large  value  of  kp/k0.  The  reduction  in  reflection  coefficient  is  less  promi¬ 
nent  near  the  critical  angle.  This  property  is  very  important  so  far  as  the  field  of  a 
dipole  is  concerned  since  a  significant  portion  of  the  contribution  of  the  integrand  to 
the  integral  comes  from  this  point.  Figure  10  shows  the  magnitude  of  the  reflection 
coefficient  for  both  polarizations  versus  normalized  surface  rms  height  for  a  canopy 
with  £i  =  1.03  +  j’0.015  at  30MHz.  Figure  11  compares  the  field  of  a  vertical  dipole 
( h  =  0.2A,  h'  =  0.1A)  obtained  using  exact  solution  and  the  approximate  distorted 
Born  solution  in  a  canopy  with  £i  =  1.01  +  j0.006  and  two  different  values,  m  =  A/4 
and  m  =  A/2  at  30MHz  where  an  excellent  agreement  is  shown.  A  similar  compari¬ 
son  is  shown  in  Fig.12  where  the  effect  of  imaginary  part  of  the  dielectric  constant  is 
examined.  The  effect  of  surface  roughness  is  shown  in  Figures  13  and  14.  Vertical  and 
horizontal(/  =  y)  dipole  in  a  canopy  with  £\  =  1.01  +  jO.006  and  £i  =  1.03  +  j0.015 
and  three  different  values  of  surface  rms  height,  a  =  0,  a  =  A/2,  and  a  =  A  are 
considered  (observation  point  at  <f>  =  45°).  It  is  shown  that  the  surface  roughness 
reduces  the  field  intensity  which  is  a  function  of  the  surface  rms  height.  It  is  very 
interesting  to  note  that  despite  relatively  significant  surface  rms  height,  lateral  wave 
is  still  the  dominant  source  of  the  field  in  a  forested  area.  This  can  be  attributed 
to  the  fact  that  most  contribution  of  the  mean  field  comes  from  the  integrand  of  the 
equation(15)  for  values  of  kp/k0  near  the  critical  angle.  As  shown  before  the  value 
of  the  reflection  coefficient  near  the  critical  angle  does  not  experience  a  significant 
reduction  due  to  the  surface  rms  height. 

5  Conclusion 

The  effect  of  canopy-air  interface  roughness  on  the  propagation  of  electomagnetic 
waves  in  forested  environment  was  investigated  in  this  paper.  An  analytic  formula¬ 
tions  that  takes  effect  of  the  rough  boundary  into  account,  are  obtained  for  both  cases 
of  plane  wave  and  an  arbitrary  oriented  dipole  excitation.  The  solution  is  obtained 
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using  distorted  Born  approximation  to  a  volumetric  integral  equation  for  the  induced 
polarization  current  in  a  hypothetical  layer  above  the  canopy.  This  formulation  was 
validated  by  comparing  the  approximate  results  with  the  exact  results  in  the  special 
case  of  smooth  interface.  It  is  shown  that  the  canopy-air  interface  roughness  reduces 
the  mean  field  surface  reflectivity  drastically  for  plane  wave  illumination  at  incident 
angles  below  the  critical  angle.  A  significant  result  of  plane  wave  simulations  was 
the  discovery  of  the  fact  that  the  mean  surface  reflectivity  near  the  critical  angle  is 
not  drastically  affected  by  the  surface  roughness,  thereby  allowing  the  propagation  of 
the  lateral  wave  despite  significant  dielectric  interface  roughness.  Direct  simulations 
of  the  field  of  an  arbitrary  dipole  in  a  forest  with  different  effective  permittivity  and 
surface  roughness  show  that  the  field  at  the  observation  point  is  anywhere  between  0 
to  lOdB  lower  than  that  for  a  smooth  air-canopy  interface  dependent  on  the  value  of 
the  rms  height  surface  roughness (&0<7). 
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Foliage 

Envelope  e2  is  replaced  with  ej 

plus  a  polarization  current 


P 


Figure  1:  A  short  dipole  embedded  in  a  forest  canopy  with  rough  interface  modeled 

by  an  effective  dielectric  constant  e\.  A  layer  above  the  canopy  is,  equivalently, 

replaced  with  a  dielectric  slab  with  planar  interface(z  =  0  plane)  and  permittivity  £i 

— *  — + 

in  addition  to  a  polarization  current,  J  =  ik\Yi{ev  —  ei)E*. 
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Distancefm] 

(b)  /  =  y 

Figure  2:  The  diffracted  field  intensity  components  (see  (15))for  a  vertical(a)  and 
horizontal(b)  dipole  in  a  canopy  with  Z\  =  1.03  +  J0.006  with  a  smooth  boundary  (a  = 
0)  ,  h  =  3m,  h'  =  3m  at  30  MHz.  The  results  are  obtained  with  the  approximation 
solution  obtained  from  the  asymptotic  integral  evaluation. 
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Figure  3:  Electric  field  of  a  vertical  dipole  as  a  function  of  distance  in  a  canopy  with 
effective  permittivity  e%.  Three  different  values  for  £\  are  used,  1.01  +  j0.006,  1.03 
+  j 0.006,  1.05  +  J0.006,  and  the  location  of  the  dipole  and  observation  points  are 
respectively,  h  =  2m  (0.2A),  h!  =  lm  (0.1A). 
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Figure  4:  Electric  field  of  a  vertical  dipole  as  a  function  of  distance  in  a  canopy 
with  effective  permittivity  ex.  Three  different  values  for  e\  are  used,  1.03  +  j0.006, 
1.03  +  j0.03,  1.03  +  j0.06,  and  the  location  of  the  dipole  and  observation  points  are 
respectively,  h  —  2m  (0.2A),  h'  =  lm  (0.1A). 


Figure  5:  Electric  field  of  a  vertical  dipole(f=30MHz)  in  a  canopy  with  effective 
dielectric  constant  e\  as  function  of  k0(h  +  h')  at  a  fixed  observation  point(p  =  2.8km). 
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(b) 

Figure  6:  Comparison  of  the  magnitude(a)  and  phase(b)  of  exact  reflection 
coefficient(i?/le2j/'um)  with  the  approximate  distorted  Born  solution(see  (8))  for  per¬ 
pendicular  polarization.  Two  permittivity  values  and  two  mean  layer  thickness  values 
are  shown  as  a  function  of  sin#,  =  kp  /  k0.  £i  —  1.03+j0.005,  and  m  =  6m  (0.6A) 
and  the  other  is  for  e\  =  1.01+j0.005,  and  m  =  10m  (A).:  frequency  is  30MHz. 
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(a) 


(b) 

Figure  7:  Comparison  of  the  magnitude(a)  and  phase(b)  of  exact  reflection 
coefficient(i?ve2jfcum)  with  the  approximate  distorted  Born  solution(see  (9))  for  par¬ 
allel  polarization.  Two  permittivity  values  and  two  mean  layer  thickness  values  are 
shown  as  a  function  of  s'm  Oi  =  kp  /  k0.  £i  =  1.03+j0.005,  and  m  =  6m  (0.6A)  and 
the  other  is  for  £i  =  1.01+j0.005,  and  m  =  10m  (A).:  frequency  is  30MHz. 
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Figure  8:  The  region  of  validity  of  the  distorted  Born  approximation  (points  under 
the  curve)  is  shown  with  regard  to  the  reflection  coefficient  for  5%  error  criterion. 
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(a) 


(b) 

Figure  9:  The  magnitude  of  reflection  coefficient  for  perpendicular(a)  and  parallel(b) 
polarizations  as  a  function  of  kp/k0  =  sin#;  for  a  canopy  with  £X  =  1.03+j0.015 
assuming  flat  surface  (cr  =  0),  a  =  A/10  ( m  =  A/5),  and  a  —  A/5  (m  =  2A/5).  It 
is  shown  that  surface  roughness  significantly  decrease  the  magnitude  of  the  reflection 
coefficient 
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Figure  10:  Magnitude  of  reflection  coefficient  versus  normalized  rms  height(&0cr)  at 
the  critical  angle  for  a  canopy  with  ex  =  1.03  +  j0.015,  and  ex  =  1.01  +  j0.005  for 
perpendicular  and  parallel  polarization. 
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Figure  11:  Field  of  a  vertical  dipole  in  a  canopy  with  ex  =  1.01+j0.006  and  smooth 
interface  obtained  from  the  exact  solution  is  compared  with  distorted  Born  approx¬ 
imation  for  two  values  of  m.  An  operation  frequency,  30MHz,  and  h  =  0.2A  and 
h'  =  0.1  A  are  assumed. 
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Figure  12:  Field  of  a  vertical  dipole  in  a  canopy  with  e\  =  1.03+j0.006  or  1.03+j0.012 
and  smooth  interface  obtained  from  the  exact  solution  is  compared  with  distorted 
Born  approximation  for  fixed  value  of  k0m  =  7t/5.  An  operation  frequency,  30MHz, 
and  h  =  0.2A  and  h!  =  0.1  A  are  assumed. 
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Figure  13:  Diffracted  field  of  a  vertical  dipole  in  a  forest  canopy  ( h  =  0.2A,  h!  =  0.1A) 
for  two  different  values  of  £X  =  1.01  +  J0.006  and  £1  =  1.03  +  j0.015  and  for  three 
different  values  of  surface  rms  height,  cr  =  0,  cr  =  A/2,  and  a  =  A. 
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Figure  14:  Diffracted  field  of  a  horizontal  dipole(/  =  y)  in  a  forest  canopy  (h  =  0.2A, 
h'  =  0.1A)  for  two  different  values  of  ei  =  1.01  +  j0.006  and  =  1.03  +  j0.015  and 
for  three  different  values  of  surface  rms  height,  a  =  0,  a  =  A/2,  and  a  =  A. 
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Abstract 

The  problem  of  wave  propagation  in  forest  is  revisited ,  In  particular ,  the  effect  of  the  non-planar  interface 
between  the  air  and  the  canopy  on  lateral  waves  is  examined .  An  analytical  formulation  is  obtained  for  the 
mean  field  when  both  the  transmitter  and  receiver  are  within  the  foliage.  This  formulation  is  based  on  distorted 
Born  approximation  and  is  shown  that  compared  to  a  planar  interface f  the  field  of  a  dipole  in  a  canopy  with 
rough  interface  is  significantly  reduced . 

1  Introduction 

A  generally  adopted  model  of  a  forest  at  HF  through  VHF,  that  attributes  wave  propagation  in  forest  to  a 
lateral  wave,  was  first  developed  by  Tamir  [1].  In  this  approach,  the  forest  is  modeled  by  a  homogeneous  half¬ 
space  dielectric  medium  with  a  planar  interface.  The  field  of  a  dipole  within  this  medium  was  then  evaluated  at 
an  observation  point  near  the  interface  using  the  asymptotic  form  of  the  integral  involved.  This  solution  shows 
that  the  field  at  the  observation  point  is  dominated  by  the  so  called  lateral  wave  that  travels  along  the  flat 
canopy-air  interface.  In  reality,  however,  the  interface  between  forest  canopies  and  air  is  not  flat,  hence  it  is  not 
clear  as  to  what  happens  to  the  lateral  wave  or  whether  it  can  even  be  excited  or  not.  In  this  paper,  the  effect 
of  roughness  of  interface  between  canopy  and  air  on  the  wave  propagation  in  forest  areas  is  investigated.  An 
analytical  solution  is  obtained  using  the  volumetric  integral  equation  in  conjunction  with  the  distorted  Born 
approximation. 

2  Analytical  Formulation 

Geometry  of  the  diffraction  problem  is  shown  in  Fig.l  where  the  dipole  located  at  heights  h  and  h 9  inside 
a  canopy  with  effective  dielectric  constant  £\.  The  envelope  of  the  canopy-air  interface  is  denoted  by  d(x,y). 
The  permittivity  of  the  upper  medium  (air)  is  denoted  by  £2.  We  modify  the  problem  by  extending  the  canopy 
to  z  =  0  plane  and  assume  that  there  exists  a  volumetric  polarization  current  J  =  ikiYifa  —  £i)E*,  where 
E*  =  Ez  +  Er  +  Es.  Here  El  is  the  incident  field,  Er  is  the  reflected  wave  from  the  planar  interface  and  E5  is 
the  scattered  field  generated  by  J  itself.  To  the  first  order  in  (s2  —  £i)  the  Born  approximation  can  be  used  to 
the  scattered  field 

_  roo  roo  rO  _ .  _ 

Es  =  fc02(e 2  -  £i)  /  /  /  G(r,  r ')  •  [E^f*)  +  Er(f')] dx'dy'dz'  (1) 

J — oo  J  —oo  J —d(x,y) 

where G  is  the  dyadic  Green’s  function  of  the  half  space  dielectric  medium  [2,  3]. 

In  (1)  d(x ,  y)  is  a  two  dimensional  random  process  describing  the  interface  between  the  canopy  and  air  and  is 
assumed  to  be  Gaussian  with  a  mean  value  of  m(a  positive  number)  and  standard  deviation  of  a.  Distorted 
Born  approximation  provides  a  more  accurate  solution  for  Es.  In  this  approximation,  a  phase  correction  term  is 
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incorporated  into  the  expressions  for  the  interna]  field  to  account  for  the  difference  in  the  propagation  constant 
between  the  air  and  the  canopy  as  seen  in  the  following  expression  of  incident  and  reflected  waves. 

E*  =  p(k\ 

(2) 

Er  =  P(— JbJ  )Rettk*x'+kvv'-k'*z'l .  c,,fc2*“fcu)TO  ; 


where  P(k\z)  is  a  unit  vector  (p  =  v  or  h)  as  is  defined  in  [2,  3]  and  RVth  is  the  Fresnel  reflection  coefficient 
for  vertical  or  horizontal  polarization. 

Since  the  layer  thickness  is  not  uniform,  the  phase  correction  terms  in  (2)  are  chosen  for  a  layer  with  a  uniform 
thickness  m.  Substituting  2-D  Fourier  transform  ofG  in  (1)  and  after  performing  the  integration  with  respect 
to  z\  evaluation  of  (E5)  would  require  computation  of  the  term  like  (e±J'srf)  which  for  a  Gaussian  d  are  found 
to  be 


•  j  esm 

W,d)  =  Vle' 


.<72S. 


,a2s* 


2  erfc{—j—j=)  +  e  2  erfc(j——)} 


<e" 


,a2s* 


■<r2s. 


-[« rH-erfci-j- +  e-^er/cy—)] 


(3) 


where  s  can  be  either  sj  =  k\z  +  k\z  or  S2  =  k2z  —  k\z. 

The  integration  with  respect  to  x'  and  y'  can  now  be  carried  out  analytically  which  results  in  S(kx-kix)-S(ky—kiy). 
This  in  turn  simplifies  the  integration  with  respect  to  kx  and  ky.  Thus  the  final  result  is  readily  obtained  as, 

E4  =  Er  +  <: Es >  «  (Rre/  +  R%lnylki**+kUv-*i*’ ]  (4) 

where  Rrej  is  the  Fresnel  reflection  coefficient  for  the  canopy-air  boundary  at  z  =  0. 

Rref  +  Rglrn  is  accurate  enough  for  horizontal  polarization,  but  it  cannot  accurately  predict  the  Brewster  angle 
for  vertical  polarization.  To  rectify  this  deficiency,  higher  order  solutions  must  be  obtained  but  it  is  sufficient 
to  use  only  term  of  the  dyadic  Green’s  function.  Thus  the  partial  second  order  solution  can  be 

obtained  and  is  given  by 


E1 


i(l.5)  u  2  (g 2  ~  gl)2 


-k0z— — — [&  •  [zz  ■  (E*  +  E r)]dv 
£2  Jv 


(5) 


The  ensemble  average  of  Es(1,5^  can  be  obtained  in  a  manner  similar  to  what  was  used  in  computation  of  (E^1)). 
After  some  algebraic  manipulations,  the  reflection  coefficients  for  the  mean  field  are  obtained  and  given  by 
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+  —  ((eis2d)  -  (e-is2d))}eis2Tn 
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Born(y) 


+  — (1  -  Rv2  +  Rv(eisid )  -  ( e~isid ))  •  (^k2  -  k[2)}eis2m 

S\  £2 


where  kp2  =  k{J  +  ky2. 

Using  the  above  result,  the  solution  for  an  infinitesimal  dipole  embedded  within  the  foliage  can  be  obtained 
by  expanding  the  field  of  the  dipole  in  terms  of  a  continuous  spectrum  of  plane  waves.  Assuming  that  a  dipole, 
whose  orientation  is  denoted  by  a  unit  vector  Z,  is  located  at  ro  =  —  (h  +  m)i,  and  using  superposition,  the 
mean  scattered  field  can  be  computed  from  the  coherent  sum  of  all  reflected  plane  waves.  That  is 


<eb„„>  =  i-[((  •  + (i  • 
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For  more  accurate  computation,  the  mean  field  for  a  dielectric  medium  with  a  rough  surface  is  calculated  by, 

Eexaci(^)  ^  Eexacf(^'  =  0)  +  E Borni E/gom^Cr  —  0) 

where  Eexaci(o  =  0)  is  the  field  of  the  dipole  in  the  presence  of  the  upper  free-space  medium  with  a  planar 
interface  at  z  =  —  m. 

3  Numerical  Simulation 


In  order  to  verify  the  accuracy  of  the  distorted  Born  approximation,  first,  the  approximate  analytical  solution 
for  reflection  coefficient  of  a  planar  boundary  is  compared  with  the  exact  Fresnel  reflection  coefficient  when  a 
plane  wave  is  incident  at  the  boundary.  We  consider  a  canopy  with  effective  permittivity  e  =  1.03+  tO.OOl 
and  the  interface  between  the  canopy  and  air  is  assumed  planar  at  a  distance  m  =  6[m]  from  the  x-y  plane 
of  the  reference  coordinate  system.  Figure  2  and  3  show  the  comparison  between  the  reflection  coefficients 
as  predicted  by  the  distorted  Born  approximation  and  the  exact  solution  for  both  parallel  and  perpendicular 
polarizations.  A  very  good  agreement  is  obtained  for  this  example.  Further  sensitivity  analysis  show  that  the 
accuracy  of  the  distorted  Born  approximation  degrades.  Similar  behavior  is  obtained  when  m  is  kept  fixed  and 
dielectric  constant  of  the  dielectric  layer  is  increased.  Next  we  considered  the  field  of  dipole  inside  a  dielectric 
layer  with  e  =  1.03  +  t'0.001  at  a  depth  of  2[m]  below  the  interface.  The  field  is  observed  at  a  depth  in  l[m] 
below  the  interface  as  a  function  of  radial  distance.  The  exact  solution  is  compared  with  the  distorted  Born 
approximation  for  a  chosen  value  of  m  =  2[m]  at  30[MHz]  in  Fig.4.  An  excellent  agreement  is  obtained.  Close 
examination  of  these  results  indicates  a  maximum  relative  error  of  2  %  between  the  two  solutions.  As  for 
the  plane  wave  illumination  simulations,  the  discrepancy  between  the  distorted  Born  and  the  exact  solution 
increases  with  increasing  m. 

With  confidence  in  the  distorted  Born  solution,  the  effect  of  the  interface  roughness  on  the  field  can  now 
be  examined.  Fig. 5  shows  the  variation  of  the  field  as  function  of  radial  distance  between  the  transmitter  and 
the  receiver  ( h  =  2[m],  h!  —  l[m])  for  two  cases.  In  the  first  case,  ey  =  1.01  +  i0.6  and  kcr  =  3  (roughness 
parameter)  and  the  second  case,  e\  =  1.03  +  «0.6  and  ka  =  2.  The  results  are  also  compared  with  those  had 
kcr  =  0  (flat  interface).  It  is  shown  that  these  surface  roughnesses  reduce  the  field  by  a  factor  of  3-5dB.  For 
these  simulations,  we  used  a  p.d.f.  for  d  of  the  following  form  fd(d)  =  y=^e_rf2/2<T2  where  d  can  only  assume 
positive  numbers. 

4  Conclusions 

Analytical  formulation  for  the  mean-field  of  a  short  dipole  embedded  in  a  forest  is  computed.  In  this 
formulation,  the  effect  of  the  roughness  of  the  air-canopy  interface  is  taken  into  account.  Distorted  Born 
approximation  is  shown  to  provide  a  very  accurate  results  for  the  limiting  case  when  the  interface  roughness 
disappears.  Simulated  results  indicate  that  the  roughness  of  the  interface  reduces  the  contribution  of  the  lateral 
waves  significantly. 
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Magnitude 


Figure  1.  Geometry  of  a  dipole  with  a  canopy 
with  rough  interface 


Phase 


Figure  2.  Comparison  of  the  magnitude  of  re¬ 
flection  coefficients  for  a  planar  interface  using 
exact  and  Born  approximation 


Figure  3.  Comparison  of  the  phase  of  reflection  Figure  4.  Comparison  of  exact  versus  approx- 

coefficients  for  a  planar  interface  using  exact  and  imate  Born  solution  for  the  field  in  a  dielectric 

Born  approximation  with  flat  interface 
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Abstract 

In  this  paper  a  preliminary  study  is  carried  out  to  demonstrate  the  application 
of  wavelets  for  improving  the  computation  time  and  reducing  computational  memory 
required  for  evaluating  the  statistics  of  the  scattered  field  from  rough  surfaces  using  the 
method  of  moments  in  conjunction  with  a  Monte  Carlo  simulation.  In  specific,  Haar  and 
the  first-order  B-spline  wavelet  basis  functions  are  applied  to  the  MoM  formulation  of 
one-dimensional  rough  surfaces  in  order  to  compare  the  computation  time  and  sparsity 
for  wavelets  in  the  same  family  but  of  higher  order.  Since  the  scattering  coefficient  (the 
second  moment  of  the  backscatter  field  per  unit  area)  is  a  gentle  function  of  the  surface 
parameters  and  the  radar  attributes,  it  is  demonstrated  that  a  relatively  high  thresholding 
level  can  be  applied  to  the  impedance  matrix  which  leads  to  a  sparser  impedance  matrix 
and  faster  computation  time.  It  is  also  shown  that  applying  a  high  threshold  level  the 
coefficients  of  the  high  order  wavelets  would  increase  out  of  proportion,  however  the  effect 
of  these  current  components  averages  out  when  computing  the  scattering  coefficients. 

The  resulting  sparse  impedance  matrices  are  solved  efficiently  using  fast  search 
routines  such  as  the  conjugate  gradient  method.  A  systematic  study  is  carried  out  to 
investigate  the  effect  of  different  threshold  levels  on  the  accuracy  versus  computing  speed 
criterion.  The  computed  scattering  coefficients  are  compared  to  previous  results  com¬ 
puted  using  a  conventional  pulse  basis  function  as  well  as  the  existing  theoretical  solutions 


1 


for  rough  surfaces.  It  is  shown  that  wavelet  basis  functions  provide  substant  ial  reductions 
in  both  memory  requirements  and  computation  time. 

I  Introduction 

The  problem  of  electromagnetic  scattering  from  rough  surfaces  has  been  the  sub¬ 
ject  of  intensive  investigation  over  the  past  several  decades  for  its  application  in  a  number 
of  important  remote  sensing  problems.  Radar  remote  sensing  of  the  oceans,  soil  mois¬ 
ture,  and  mine  detection  using  wideband  radars  are  such  examples.  For  these  problems, 
where  the  rough  surface  is  either  the  primary  target  or  the  clutter,  the  understanding 
of  interaction  of  electromagnetic  waves  with  the  rough  surface  is  essential  for  developing 
inversion  or  detection  algorithms.  An  exact  analytical  solution  for  random  rough  surfaces 
does  not  exist.  However,  approximate  analytical  solutions  exists  for  rough  surfaces  with 
specific  types  of  surface  roughness  conditions.  For  surfaces  with  small  root  mean  square 
(rms)  height  and  slope,  the  small  perturbation  method  (SPM)  is  the  most  commonly 
used  formalism.  Formulations  based  on  SPM  exist  for  perfectly  conducting  [10],  homoge¬ 
neous  dielectric  [4],  and  inhomogeneous  dielectric  [12]  rough  surfaces.  Another  classical 
solution  which  is  valid  for  surfaces  with  large  radii  of  curvature  is  based  on  the  tangent 
plane  approximation  [15].  The  region  of  validity  of  these  classical  approaches  are  rather 
limited.  In  recent  years  much  effort  has  been  devoted  to  extend  the  region  of  validity  of 
these  models  [3,  9],  however,  the  improved  techniques  still  have  the  basic  limitations  of 
the  original  models. 

An  alternative  approach  for  evaluating  of  the  scattered  field  and  its  statistics  for 
rough  surfaces  is  Monte  Carlo  simulation.  In  this  approach  many  sample  surfaces  with 
the  desired  roughness  statistics  are  generated,  and  then  the  scattering  solution  for  each 
sample  surface  is  obtained  using  a  numerical  method.  Monte  Carlo  simulation  have 
primarily  been  considered  for  evaluating  performance  of  and  characterizing  the  region  of 
validity  of  approximate  analytical  models  [1,  2,  3,  9].  In  general,  the  limitations  of  Monte 
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Carlo  simulation  of  scattering  from  rough  surfaces  are  the  computation  time  and  the 
required  memory  as  the  typical  size  of  the  scatterers  (sample  surfaces)  must  be  chosen 
to  be  much  larger  than  the  wavelength.  Another  issue  is  that  the  rough  surfaces  are 
the  targets  of  infinite  extent  which  must  be  truncated  appropriately  before  the  numerical 
scattering  solution  can  be  obtained.  This  can  be  done  either  using  a  tapered  illumination 
[9],  or  padding  the  sample  surfaces  with  a  tapered  resistive  sheet  [11].  It  has  been  shown 
that  with  the  tapered  illumination,  larger  sample  surfaces  must  be  used  as  a  considerable 
portion  of  the  induced  currents  on  the  surface  do  not  contribute  significantly  to  the 
total  scattered  field.  Application  of  the  tapered  resistive  sheet  is  advantageous  in  that 
a  relatively  small  portion  of  the  sample  surfaces  is  actually  used  to  suppress  the  edge 
currents.  This  improves  the  computation  time  and  reduces  the  required  memory. 

In  order  to  use  Monte  Carlo  simulation  for  evaluating  the  scattering  statistics  of 
rough  surfaces  more  routinely,  computationally  more  efficient  scattering  codes  must  be 
developed.  In  this  paper  the  application  of  wavelets  as  a  basis  function  for  the  expan¬ 
sion  of  induced  surface  currents  is  considered.  Traditional  method  of  moments  (MoM) 
in  conjunction  with  Galerkin’s  method  would  require  matrix  fill  computation  time  of 
the  order  of  N 2  and  matrix  inversion  computation  time  of  the  order  of  N3  (using  Gaus¬ 
sian  elimination).  It  is  well  known  that  the  solution  of  linear  system  of  equations  can 
be  obtained  far  more  efficiently  using  search  routines,  such  as  the  Conjugate  Gradient 
method,  if  the  matrix  of  the  coefficients  is  a  sparse  matrix.  In  MoM,  the  application 
of  conventional  pulse  or  rooftop  basis  and  testing  functions  would  usually  produce  full 
impedance  matrices.  Although  the  diagonal  elements  are  usually  larger  than  the  rest  of 
the  elements,  the  smaller  elements  cannot  be  arbitrarily  thresholded  without  drastically 
altering  the  resulting  scattering  pattern.  The  success  of  wavelet  expansion  function  in 
generating  sparse  matrices  have  been  demonstrated  for  many  circuits  and  antenna  prob¬ 
lems  [5,  6,  7,  8].  In  the  Monte  Carlo  simulation  of  scattering  from  rough  surfaces  the 
quantities  of  interest  are  the  statistical  parameters,  such  as  the  mean  and  variance  of  the 
scattered  field,  and  therefore  it  is  expected  that  the  overall  accuracy  be  less  sensitive  to 
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the  threshold  level. 


An  investigation  is  conducted  on  the  use  of  two  different  types  of  wavelets  with 
compact  support,  Haar  and  B-spline  wavelets  with  edge  wavelets,  and  the  effect  of  dif¬ 
ferent  threshold  level  with  regard  to  the  overall  accuracy  and  the  computation  time. 
The  method  is  applied  to  one-dimensional  perfectly  conducting  random  rough  surfaces 
to  demonstrate  the  improvements  achieved.  In  the  Monte  Carlo  analysis  presented  here 
the  tapered  resistive  sheet  approach  is  used  to  suppress  the  edge  current  for  plane  wave 
illumination.  The  numerical  results  are  also  compared  with  the  approximate  analytical 
solutions. 


II  Integral  Formulation  of  Scattering  Problem 

In  order  to  characterize  both  the  backscattering  coefficient  and  the  bistatic  scatter¬ 
ing  coefficient  using  the  MoM  and  a  Monte  Carlo  simulation,  a  large  number  of  random 
surfaces  with  known  statistical  parameters  ( ks ,  M)  must  be  generated.  Then  it  is  desired 
to  find  the  surface  current  density  Je  induced  by  a  plane  wave  from  which  the  scattered 
field  can  be  computed.  For  a  horizontally  polarized  plane  wave  excitation  the  scattered 
electric  field  is  given  by: 

E’(p)  = JtMi>)H!,'\k0\p -?\)de  ,  (i) 

here  ko  is  the  wave  number,  Zq  is  the  intrinsic  impedance  of  free  space,  is  the  zeroth- 

order  Hankel  function  of  the  first  kind,  and  p  and  pf  are  the  position  vectors  of  observation 
and  source  points,  respectively.  The  sample  surface  is  discretized  into  sufficiently  small 
cells,  as  shown  in  Figure  1,  and  equation  (1)  is  cast  into  a  matrix  equation.  The  expression 
for  the  plane  wave  propagating  along  k{  =  sin  0{X  —  cos  is  given  by, 


E  (#m,  Dm) 


eik0(xm  sin  Oi  —ym  cos  0») 


(2) 
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As  mentioned  before,  because  of  the  singular  behavior  of  the  current  near  the  edges 
of  the  surface  when  excited  by  a  horizontally  polarized  incidence  wave,  tapered  resistive 
sheets  must  be  added  at  the  edges  of  the  sample  surface  in  order  to  suppress  the  edge 
currents.  The  induced  surface  current  on  a  resistive  sheet  is  proportional  to  the  tangential 
electric  field,  or  mathematically: 


fix  (fix  E)  =  -R3  .  (3) 

where  R  is  the  surface  resistivity  (for  a  perfect  conductor,  R  =  0).  Another  boundary 
condition  for  resistive  sheets  mandates  continuity  of  the  tangential  electric  field  across 
the  resistive  sheet,  that  is,  [n  x  E]+  =  0. 

Therefore  the  electric  field  integral  equation  for  the  surface  current  is  given  by 

K(p)  =  mMM )  +  ^  j 3 *W'))Hl'\h,\m  -  7(C) \)dt  .  (4) 

III  Tapered  Resistive  Sheets 

Tapered  resistive  sheets,  as  introduced  in  section  II,  are  added  on  the  edges  of  the 
surface  samples  to  suppress  edge  effects  on  the  induced  surface  current  and  the  scattered 
fields.  An  optimum  tapered  function  for  resistivity  (found  by  trial  and  error)  is  given 
by  [11]: 


o,  4  |*|  <  § 


(5) 


where  D  is  the  width  of  the  sample  surface  and  Dp>  is  the  width  of  the  resistive  section. 
This  taper  function  was  found  to  significantly  decrease  diffraction  due  to  the  edge  discon¬ 
tinuity.  Figure  2  shows  a  normalized  resistivity  profile  of  the  tapered  resistive  sheets  and 
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the  placement  of  the  resistive  sheets  on  a  surface  described  by  a  Gaussian  hump.  Fig¬ 
ures  3a  and  3b  show  the  induced  current  distribution  on  a  flat  surface  with  and  without 
thin  tapered  resistive  sheets,  respectively. 

Expanding  the  current  in  terms  of  the  basis  functions  Je{p)  =  Y.an<f>n{p)  and  ap¬ 
plying  Galerkin’s  method  to  equation  (4)  we  have: 


£  ^  jJanK(i>)K(p')H$\h\p  -  ?\ )dfd(\  ,(6) 

where  <j>m(p )  is  the  testing  function,  and  <t>n{p')  is  the  basis  function.  Equation  (6) 
can  be  solved  using  numerical  integration,  and  cast  into  a  matrix  equation,  given  by 
[Z]  [I]  =  [V],  This  matrix  equation  can  easily  be  solved  to  find  the  surface  current  Je. 

IV  Applications  of  Compact-Support  Wavelets 

As  stated  above,  the  EFIE  used  in  conjunction  with  Galerkin’s  method  can  be  cast 
into  a  matrix  equation  using  appropriate  expansion  and  testing  functions.  The  restric¬ 
tion  on  the  expansion  and  testing  functions  is  that  they  have  to  be  in  the  domain  of  the 
integral  equation  operation.  To  expand  the  induced  current  in  terms  of  a  multiresolu¬ 
tion  expansion,  first,  the  current  must  be  projected  onto  the  x-axis  or  the  surface  must 
be  arclength  parameterized.  For  natural  rough  surfaces  with  moderate  rms  slope  it  is 
more  convenient  to  project  the  current  on  the  x-axis  since  the  domain  will  be  identical 
for  all  sample  surfaces.  Applying  a  multiresolution  expansion  to  the  projected  current 
distribution  ( /(x ))  on  the  x-axis  we  have: 

mh-l 

/(x)  =  y  =  ^  “I"  ^  ^  (7) 

n  n  m=mi  n 
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which  consists  of  the  scaling  functions  at  the  lowest  resolution,  plus  wavelets  at  the 
lowest  resolution  and  subsequent  higher  resolutions.  This  expansion  is  equivalent  to  one 
consisting  of  only  scaling  functions  at  the  highest  resolution  (<f>mhiTl(x)).  Using  the  Fast 
Wavelet  Transform(FWT)[18],  the  program  needs  only  to  compute  the  EFIE  impedance 
matrix  at  the  highest  level  of  resolution.  From  this  computed  impedance  matrix  at  the 
highest  resolution,  lower  resolutions  can  be  found  in  terms  of  the  original  computed 
impedance  matrix.  This  reduces  computation  time  drastically  from  having  to  perform 
the  integration  for  each  wavelet  at  various  resolution  levels. 

Sparse  matrices  arise  due  to  the  fact  that  the  scaling  functions  and  wavelets  are 
orthogonal  and  wavelets  have  zero  moments.  For  example,  Haar  and  linear  B-spline 
wavelets  both  have  zero  mean  and  the  first  moment  of  the  linear  B-spline  wavelet  is  also 
zero.  Higher  order  wavelets  also  have  vanishing  higher  moments,  that  is 

J  xnip(x)  dx  =  0  n  €  {0,...,N},  (8) 

where  N  depends  on  the  order  of  the  B-spline  wavelet.  Figures  4  and  5  show  the  Haar 
and  linear  B-spline  scaling  functions  and  wavelets. 

The  significance  of  vanishing  moments  of  the  wavelets  stems  from  the  fact  that 
the  integration  of  the  kernel  over  the  domain  of  the  wavelet  would  produce  a  small 
quantity.  This  quantity  is  smaller  for  wavelets  with  higher  vanishing  moments.  Then, 
when  integrating  over  the  kernel,  very  small  matrix  element  values  are  usually  calculated 
for  cells  that  are  relatively  far  from  one  another  with  the  kernel  being  fairly  regular. 
For  adjacent  cells  and  self-cells  the  integration  of  the  kernel  produce  element  values  that 
are  significant.  In  fact,  the  self-cell  contribution  is  the  largest  element  of  the  impedance 
matrix. 

Once  the  impedance  matrix  for  the  highest  resolution  is  computed,  the  FWT  algo¬ 
rithm  is  applied  to  the  matrix  to  find  an  equivalent  impedance  matrix  for  the  multireso¬ 
lution  expansion.  Using  the  FWT  it  is  expected  that  a  sparse  matrix  will  be  generated. 


7 


At  this  stage  a  user-defined  threshold  level,  usually  of  the  order  of  .01%  to  1%,  is  imposed 
on  the  matrix  and  only  significant  elements  of  the  matrix  will  be  preserved. 

For  a  preliminary  examination  of  the  wavelet-based  Method  of  Moment,  a  Gaussian 
hump  described  by  equation 


y(x )  =  A  •  e  (  u*° )  (9) 

is  used.  The  effects  of  the  multiresolution  expansion  and  application  of  different  threshold 
levels  on  the  moment  matrix  on  the  bistatic  scattering  is  investigated.  Initially,  pulse 
basis  functions  with  the  highest  resolution  is  used  to  generate  a  reference  solution  to  the 
problem  of  a  plane  wave  incident  at  normal  incidence  upon  a  Gaussian  hump  surface 
and  the  solution  is  found  with  no  threshold  applied.  This  solution  became  the  standard 
of  comparison  for  subsequent  tests  on  the  Gaussian  hump  surface  involving  both  Haar 
wavelets  and  first  order  B-spline  wavelets  for  various  levels  of  resolution  and  threshold 
levels. 

Using  Haar  wavelets,  described  by  equation 


|  1,  0  <  |z|  <  \ 

\  “!  \  <  1*1  <  1 

along  with  pulse  basis  functions,  the  moment  matrix  produced  had  a  maximum  sparsity 
level  of  over  89%  when  a  threshold  level  of  0.3%  was  applied  and  5  levels  of  resolution 
were  used.  The  bistatic  scattering  pattern  produced  from  the  reduced  moment  matrix 
shows  no  significant  differences  in  the  scattering  levels  when  compared  with  the  results 
obtained  from  the  original  full  moment  matrix.  The  sparsity  level  increases  with  the 
number  of  resolution  levels,  because  a  larger  number  of  localized  wavelets  contribute 
to  the  cancellation  effect  at  slowly  varying  regions  of  the  induced  current  distribution. 
Sharply  varying  components  of  the  current  which  are  mostly  localized  are  captured  by  a 
small  number  of  wavelets.  Keeping  the  same  threshold  level  for  the  moment  matrix,  the 
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sparsity  achieved  for  2,  3,  and  4  levels  of  resolution  are  respectively  72.5%,  80.7%).  and 
84.8%. 

Next,  rooftop  basis  functions  along  with  linear  B-spline  wavelets  w^ere  investigated. 
The  results  from  the  linear  B-spline  wavelet  case  agreed  well  wdth  the  Haar  case,  and  a 
sparsity  of  over  99%  was  observed  w'hen  a  threshold  level  of  3%  w'as  applied  and  5  levels 
of  resolution  were  used,  and  the  overall  bistatic  scattering  pattern  shows  no  significant 
differences.  As  expected,  higher  sparsity  is  achieved  with  linear  B-spline  wavelets.  For  2, 
3,  and  4  levels  of  resolution  with  the  same  threshold  level  applied  to  the  moment  matrix, 
sparsities  of  93.2%,  97.3%,  and  97.8%,  respectively,  are  achieved  for  the  linear  B-spline 
wavelet  expansion.  Once  a  sparse  matrix  is  obtained,  a  search  routine,  such  as  conjugate 
gradient  method,  can  be  used  to  find  the  solution.  Since  the  number  of  non-zero  matrix 
elements  is  far  less  than  the  original  matrix,  the  computation  time  for  obtaining  the 
solution  is  reduced  drastically.  The  comparison  between  the  full  matrix  solution,  Haar, 
and  linear  B-spline  wavelets  with  5  levels  of  resolution  and  the  threshold  levels  stated 
above  for  the  respective  basis  functions  is  shown  in  Figure  6. 

V  Monte  Carlo  Simulations  of  Random  Rough  Sur¬ 
faces 

Encouraged  by  results  obtained  for  a  simple  Gaussian  hump  surface,  the  multi¬ 
resolution  expansion  method  is  then  applied  to  random  rough  surfaces  of  known  statistical 
parameters.  Near  normal  incidence,  sample  surfaces  for  numerical  analysis  must  be 
at  least  40  correlation  lengths  (£)  long  in  order  to  accurately  characterize  the  bistatic 
pattern  [13].  This  requirement  becomes  more  stringent  for  near  grazing  incidence. 

Since  the  Monte  Carlo  simulation  of  the  scattering  problem  requires  the  numerical 
solution  of  the  problem  many  times,  the  computation  speed  achieved  by  thresholding 
the  moment  matrix  of  a  multi-resolution  expansion  in  conjunction  with  a  search  routine 
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linear  system  solver  becomes  very  significant  with  regard  to  the  overall  computation  time 
necessary  for  evaluating  the  statistics  of  the  scattered  field. 

V.l  Random  Surface  Generation 

Monte  Carlo  simulations  require  a  large  number  of  sample  surfaces  of  a  random 
process  with  prescribed  surface  height  statistics.  To  generate  the  sample  surfaces,  the 
procedure  in  [10,  15]  is  used.  First,  a  long  string  of  numbers  is  generated  using  a  ran¬ 
dom  number  generator  having  the  same  pdf  as  the  height  distribution  of  the  surface 
(for  example,  a  zero-mean  Gaussian  pdf).  Next,  a  subset  of  the  numbers  of  the  string 
are  correlated  with  a  weight  vector  related  to  the  Fourier  transform  of  the  desired  au¬ 
tocorrelation  function  [9].  For  the  simulations  presented  in  this  paper,  surfaces  with 
Gaussian  correlation  function  and  Gaussian  height  distribution  are  considered.  Hence, 
the  surface  statistics  are  uniquely  specified  by  the  surface  height  standard  deviation  (rms 
height),  s,  and  by  the  surface  correlation  length,  t.  Figure  7  shows  a  sample  surface  of 
a  Gaussian  process  with  ks  =  0.3,  k£  =  3.0  and  the  corresponding  histogram  of  height 
and  calculated  correlation  function  generated  from  60  independent  sample  surfaces.  The 
calculated  correlation  function  gives  kl  =  3.17,  and  the  calculated  standard  deviation 
of  height  distribution  gives  ks  =  0.3020.  These  agree  closely  with  the  desired  surface 
roughness  parameters. 

V.2  Validation  and  Results 

Numerical  simulation  of  rough  surface  scattering  is  performed  for  three  different 
surfaces  denoted  by  Si,  S2,  and  S3.  The  roughness  parameters  (ks,  kl)  for  each  of  these 
surfaces  are  respectively  ksi  =  0.3,  kl\  =  3.0;  ks2  =  0.5,  kl2  =  6.13;  and  ks3  =  2.0,  kl3  = 
2.5.  The  first  two  surfaces  fall  within  the  region  of  validity  of  small  perturbation  and 
physical  optics  models  respectively,  and  hence  the  two  numerical  results  can  be  compared 
with  the  analytical  models.  Comparisons  are  also  made  on  the  threshold  level  imposed 
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on  the  moment  matrix  for  both  Haar  and  B-spline  wavelet  basis  functions  and  on  the 
matrix  solving  times  using  a  fast  Conjugate  Gradient  solver  for  sparse  matrices.  Another 
test  that  was  run  was  the  effect  of  the  MoM  and  wavelet  technique  on  backscattering 
enhancement.  Finally,  the  effect  of  the  number  of  resolution  levels  on  the  scattering 
pattern  is  investigated. 

SPM  is  known  to  be  valid  when  ks  <  0.3,  k(  <  3.0,  and  m  <  0.3,  where  m  is  the 
rms  slope  and  is  given  by  m  —  \/2s/£  for  a  surface  with  a  Gaussian  correlation  function. 
The  analytical  bistatic  pattern  for  the  SPM  is  derived  in  [17]  for  a  3-dimensional  surface 
with  a  Gaussian  correlation  function.  For  a  2-dimensional  surface,  the  bistatic  scattering 
coefficient  (cr^)  is  related  to  the  3-dimensional  bistatic  scattering  coefficient  (<j|d)  via 
a2d  =  f  <?3d-  Then,  cr°d  is  given  by: 

°2d{0ii0s)  =  4k3 cos2(Os)cos2(0i)fhhW(\k1  -  (11) 

where  W(k)  is  the  power  spectral  density.  For  a  Gaussian  surface  correlation,  W{\k±  —  fcj_tj) 
is  given  by 


Wflfcx  -  M)  =  (12) 

Furthermore,  for  a  perfect  electrical  conductor  (PEC),  fhh  =  cos2(<f>s  —  o% )  which  equals 
1  for  a  2-dimensional  problem. 

The  incoherent  bistatic  scattering  coefficient  for  0\  =  30°  using  the  Monte  Carlo 
simulation  with  40  independent  samples  and  threshold  applied  is  compared  to  the  ana¬ 
lytical  bistatic  SPM  from  equation  (11)  in  Figure  8,  and  a  good  agreement  is  observed. 

The  Physical  optics  (PO)  method  using  the  tangent-plane  technique  for  approxi¬ 
mating  the  fields  on  a  surface,  S,  is  next  investigated.  Under  the  tangent-plane  technique, 
the  fields  present  at  any  point,  P,  on  S  are  approximated  by  the  fields  that  would  be 
present  on  a  plane  tangent  to  P.  This  is  a  valid  approximation  if  every  point  on  S  has  a 
large  radius  of  curvature.  The  three-dimensional  bistatic  scattering  coefficient  is  derived 
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in  [19]  for  the  PO  approximation.  Using  this  formulation  cr^d  *s  calculated  and  is  given 
by: 


<?u{6x,0s)  =  k 


1  -f  cos(0t)  cos(Os)  —  sin(0,)  sin(05)' 
cos(0,)  +  cos(0s) 


J 


(13) 


where 
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where  x  =  ks[cos(0i)  +  cos(0s)].  For  surface  S2  which  falls  into  the  PO  region  of  validity, 
a  comparison  is  made  in  Figure  9  between  the  PO  model  in  equation  (13)  and  results 
obtained  using  the  wavelet  based  moment  method.  The  bistatic  scattering  results  from 
S2  agree  well  with  the  theoretical  Physical  Optics  solution. 


A  comparison  of  different  threshold  levels  imposed  on  the  moment  matrix  for  surface 
Si  using  B-spline  basis  functions  is  shown  in  Figure  10.  The  parameters  for  the  simulation 
on  Si  are  0,-  =  30.0°,  length  =  32A,  resistive  tapered  ends  length  =  1A  and  Ax  =  0.1A. 
For  this  simulation,  only  a  single  surface  (N  =  1)  is  used  for  comparison.  As  is  shown 
in  Figure  10,  the  scattering  pattern  varies  very  slightly  and  only  at  angles  near  grazing 
observation.  This  figure  indicates  that  a  sparsity  of  more  than  90%  can  be  achieved 
without  substantial  compromise  in  the  accuracy  of  the  bistatic  pattern  using  B-spline 
wavelets. 


The  bistatic  scattering  coefficient  obtained  using  the  Haar  and  B-spline  wavelet  of 
surface  S3  are  shown  in  Figure  11.  For  this  simulation,  ks  =  2.0,  kt  =  2.5,  number  of 
independent  surfaces  N  =  40,  sample  length  =  32A,  resistive  tapered  ends  length  =  1A, 
Ax  =  0.04A  and  0t  =  30.0°.  Since  this  surface  does  not  fall  into  the  region  of  validity 
for  any  analytical  models,  no  comparison  may  be  made.  This  figure  shows  again  that  by 
thresholding  the  moment  matrix,  a  high  sparsity  can  be  achieved  without  compromising 
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the  accuracy.  It  is  also  shown  that  a  higher  sparsity  is  achieved  using  the  linear  B- 
spline  wavelet.  An  average  sparsity  of  97.3%  is  achieved  by  the  linear  B-spline  wavelet 
whereas  using  Haar  wavelets  only  82.2%  average  sparsity  is  obtained  for  a  threshold 
level  of  0.15%.  Figure  12  shows  a  comparison  of  the  exact  current  distribution  and  one 
where  the  impedance  matrix  has  a  sparsity  of  about  97%.  A  portion  of  a  surface  with 
characteristics  given  above  (ks  =  2.0,  k£  =  2.5,  etc.)  was  used  and  the  sparse  matrix  has 
a  threshold  of  0.5%.  Even  though  this  sample  surface  (which  is  indicitive  of  the  current 
distribution  for  every  sample  surface)  has  a  current  distribution  that  varies  significantly 
from  the  exact  solution,  because  the  far-field  scattering  pattern  is  an  averaging  process 
the  bistatic  scattering  pattern  does  not  vary  significantly  from  Figure  11. 

The  threshold  level  imposed  on  the  moment  matrix  and  the  number  of  multiplica¬ 
tions  needed  to  solve  the  matrix  using  a  fast  Conjugate  Gradient  solver  routine  is  studied 
next.  As  stated  previously,  by  imposing  a  threshold  level  the  scattering  pattern  remains 
relatively  unchanged,  yet  the  moment  matrix  could  be  made  quite  sparse.  It  was  also 
found  that  linear  B-spline  wavelets  would  produce  a  sparser  matrix  for  a  given  threshold 
level  before  the  scattering  pattern  were  to  deviate  significantly  from  the  exact  solution. 
This  is  due  to  the  fact  that  the  linear  B-spline  wavelet  has  a  vanishing  first  moment. 
Table  1  provides  the  average  number  of  multiplications  necessary  to  solve  for  the  sur¬ 
face  current  of  S3  for  both  Haar  and  linear  B-spline  wavelets  and  for  different  values  of 
threshold  levels.  It  is  found  that  the  number  of  multiplications,  or  equivalently  the  com¬ 
putation  time,  decreases  significantly  with  the  first  order  B-spline  scaling  function  and 
its  corresponding  wavelets.  For  the  Haar  wavelet-based  MoM,  a  slight  improvement  was 
observed  (approximately  20%),  yet  for  the  linear  B-spline  wavelet-based  MoM,  a  factor 
of  about  20  improvement  was  observed. 

The  question  arises  if  the  solution  produced  using  a  wavelet  based  MoM  technique 
show  the  effect  of  backscattering  enhancement.  The  backscattering  enhancement  is  pro¬ 
duces  primarily  due  to  multi-path  and  surface-wave  effects  on  a  rough  surface.  A  surface 
was  chosen  that  has  physical  parameters  where  backscattering  enhancement  is  known 
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to  occur,  and  these  parameters  are  ks  =  10.636,  k(  =  19.472.  Figure  13  was  produces 
using  the  aforementioned  physical  parameters,  as  well  as  number  of  independent  surfaces 
N  =  400,  sample  length  =  30A,  resistive  tapered  ends  length  =  1A,  Ax  =  0.0625A  and 
6i  =  10.0°.  Also,  5  levels  of  resolution  were  used,  and  the  matrices  had  an  average  spar¬ 
sity  of  97.5%.  It  can  be  seen  that  around  6S  =  —10°  there  is  an  enhanced  backscattering 
effect.  This  figure  is  in  linear  scale  to  display  the  enhanced  backscattering  effect,  and 
agrees  very  well  with  Figure  5  of  reference  [20]. 

The  effect  of  the  number  of  resolution  levels  on  the  scattering  pattern  is  next  in¬ 
vestigated  using  surface  S3.  Since  it  has  already  been  determined  that  both  Haar  and 
B-spline  wavelet-based  MoM  produce  similar  results,  the  effect  of  the  number  of  resolu¬ 
tion  levels  is  demonstrated  with  Haar  wavelets  only.  Monte  Carlo  simulations  are  run  on 
random  rough  surfaces  for  5,  4,  and  2  levels  of  resolution,  for  threshold  levels  of  0.1%  and 
0.3%,  except  for  the  Haar  with  5-levels  of  resolution  in  which  0.1%  and  0.15%  threshold 
level  is  used.  This  is  because  0.3%  threshold  level  is  too  high  for  5  levels  of  resolution  and 
the  conjugate  gradient  solver  does  not  converge.  The  results  based  on  the  full  matrix  is 
also  included  at  5  levels  of  resolution  for  comparison.  As  is  shown  in  Figure  14  the  scat¬ 
tering  pattern  starts  deviating  from  the  exact  pattern  for  5  levels  of  resolution  and  0.1% 
threshold  level  imposed  on  the  moment  matrix.  The  average  sparsity,  obtained  from  40 
independent  samples,  for  each  case  in  Figure  14  is  summarized  in  Table  2.  The  scatter¬ 
ing  pattern  agrees  quite  well  with  the  exact  solution  for  almost  all  levels  of  resolutions 
shown  with  only  slight  deviation  at  near  grazing  observations.  One  notable  exception 
is  for  4  levels  of  resolution  and  an  imposed  tolerance  level  of  0.3%,  which  deviates  no¬ 
ticeably  from  the  exact  bistatic  pattern  at  certain  angles  of  observation.  As  shown  in 
Table  2,  the  matrix  is  less  sparse  for  fewer  levels  of  resolution  at  a  single  stated  threshold 
level.  Therefore,  by  increasing  the  number  of  resolution  levels,  while  holding  a  constant 
threshold  level,  the  sparsity  of  the  matrix  will  increase,  as  is  shown  in  Table  2. 


14 


VI  Conclusion 


It  has  been  found  that  using  wavelet  basis  functions  with  MoM  and  Galerkin’s 
method  along  with  a  fast  solver  routine  such  as  conjugate  gradient  can  drastically  reduce 
both  the  memory  requirements  of  a  system  and  the  time  necessary  to  solve  the  MoM 
matrix.  This  leads  to  solutions  for  rough  surface  scattering  that  are  quite  accurate  when 
compared  to  other  basis  functions,  yet  take  significantly  less  memory  and  time  to  solve. 
Matrices  can  be  made  over  97%  sparse  yet  still  produce  accurate  bistatic  scattering 
coefficients  in  scattering  problems.  Thus,  it  becomes  possible  to  generate  statistics  for 
the  scattering  from  surfaces  of  different  roughness  in  a  relatively  short  period  of  time. 

The  number  of  resolution  levels  were  shown  to  play  a  significant  role  in  determining 
the  sparsity  of  the  matrix  and  the  accuracy  of  the  solution.  It  was  shown  that  the 
higher  the  number  of  resolutions,  the  more  sparse  the  matrix  could  be  made  without 
compromising  the  bistatic  scattering  pattern.  Also,  the  higher  order  the  B-spline  was 
made,  the  higher  the  sparsity  achieved  in  the  matrix,  and  thus,  the  faster  the  computation 
time  to  solve  the  matrix. 
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2.5 


Surface  current  -  tapered  resistive  sheet 


Figure  3:  Surface  currents  on  a  32-A  flat  plate  when  a)resistive  tapered  sheets  are  used 
and  b)  resistive  tapered  sheets  are  not  used.  The  surface  (not  shown)  is  a  PEC  flat  sheet 
with  normal  incidence. 


Figure  4:  Pulse  basis  function  and  the  corresponding  Haar  wavelet 
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Figure  6:  Comparison  between  the  exact  solution,  Haar  wavelets  with  5  levels  of  res¬ 
olution,  and  B-spline  wavelets  with  5  levels  of  resolution.  Normal  incidence  is  used  on 
a  surface  described  by  a  Gaussian  hump.  Threshold  levels  are  0.3%  and  3.0%  for  Haar 
and  B-spline  moment  matrix  with  sparsity  levels  of  89.6%  and  99.1%  respectively. 
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Figure  7:  (a)  Generated  random  surface,  with  mean  surface  height  shown  (b)  Theoret¬ 
ical  and  Monte  Carlo  probability  density  function  of  the  surface  height  distribution,  ks 
(c)  Theoretical  and  Monte  Carlo  correlation  function 


Angle  (in  Degrees) 

Figure  8:  Comparison  of  SPM  with  Monte  Carlo  simulation.  Generated  surface  param¬ 
eters  are:  ks  =  0.3,  kl  =  3.0,  N  =  40,  length  =  32A,  resistive  tapered  ends  length  =  1A, 
Ax  =  0.1A,  and  8t  =  30.0°. 
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Figure  9:  Comparison  of  the  Physical  Optics  model  and  a  random  surface  with  ks  =  0.5, 
kl  =  6.13,  N  =  50,  length  =  100A,  resistive  tapered  ends  length  =  1  A,  Ax  =  0.2A  and 
0i  =  30.0°. 


Figure  10:  Bistatic  Scattering  from  a  single  random  surface  with  ks  =  0.3,  kt  =  3.0, 
N  =  1,  length  =  32A,  resistive  tapered  ends  length  =  1A,  Ax  =  0.1A  and  6,  =  30.0°.  The 
various  threshold  levels  imposed  on  the  moment  matrix  and  its  corresponding  sparsity 
level  are  given  in  the  figure. 


23 


r 


Angle  (in  Degrees) 

Figure  11:  Comparison  of  threshold  levels  on  scattering  pattern  for  a  random  surface 
with  ks  =  2.0,  kl  —  2.5,  N  —  40,  length  =  32A,  resistive  tapered  ends  length  =  1A, 
Ax  =  0.04A,  6i  =  30.0°,  and  5  levels  of  resolution. 


Position  (in  A.) 

Figure  12:  Comparison  of  exact  current  distributions  with  the  current  distribution  from 
an  impedance  matrix  that  is  about  97%  sparse  (using  B-spline  wavelets  and  a  threshold 
of  0.5%).  A  portions  of  a  random  surface  with  ks  =  2.0,  kt  =  2.5,  length  =  32A,  resistive 
tapered  ends  length  =  1A,  Ax  =  0.04A,  —  30.0°,  and  5  levels  of  resolution  is  plotted. 


24 


Angle  (in  Degrees) 


Figure  13:  Test  for  backscatter  enhancement  effect.  The  surfaces  under  test  have 
ks  =  10.636,  kt  =  19.472,  N  =  400,  length  =  30A,  resistive  tapered  ends  length  =  1A, 
Ax  =  0.0625A  and  0,  =  10.0°.  As  can  be  seen,  a  significant  backscatter  enhancement  is 
shown  at  0S  =  —10°. 


Angle  (in  Degrees) 

Figure  14:  Comparison  of  resolution  levels  used  on  scattering  pattern  for  a  random 
surface  with  ks  =  2.0,  kl  —  2.5,  N  =  40,  length  =  32A,  resistive  tapered  ends  length  =  1A, 
Ax  =  0.04A  and  0t-  =  30.0°. 
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Threshold 

Sparsity 

Number  of  Mults 

Bspline 

0.0% 

0.0% 

5.3  x  10s 

Bspline 

0.05% 

90.7% 

8.4  x  107 

Bspline 

0.1% 

93.4% 

6.4  x  107 

Bspline 

0.5% 

97.3% 

2.6  x  107 

Haar 

0.0% 

0.0% 

7.0  x  108 

Haar 

0.05% 

66.1% 

9.5  x  108 

Haar 

0.1% 

76.9% 

7.0  x  108 

Haar 

0.15% 

82.2% 

5.6  x  108 

Table  1:  Number  of  multiplications  needed  to  solve  the  matrix  and  sparsity  vs.  imposed 
threshold  level  for  a  random  surface  with  ks  =  2.0,  kl  =  2.5,  N  =  40,  length  =  32A, 
resistive  tapered  ends  length  =  1A,  Ax  =  0.04A,  =  30.0°,  and  5  levels  of  resolution. 


Threshold 

Sparsity 

Haar 

0.0% 

0.0% 

Haar-5  levels  res. 

0.1% 

76.9% 

Haar-5  levels  res. 

0.15% 

82.2% 

Haar-4  levels  res. 

0.1% 

68.9% 

Haar-4  levels  res. 

0.3% 

85.5% 

Haar-2  levels  res. 

0.1% 

51.1% 

Haar- 2  levels  res. 

0.3% 

74.4% 

Table  2:  Sparsity  vs.  threshold  for  a  number  of  different  resolution  levels  used  for  a 
random  surface  with  ks  —  2.0,  kl  =  2.5,  N  =  40,  length  =  32A,  resistive  tapered  ends 
length  =  1A,  Ax  =  0.04A  and  0,  =  30.0°. 
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Abstract 

The  application  of  iterative  Physical  Optics  (PO)  in  conjunction  with  a  Monte  Carlo  simulation  for 
characterizing  the  bistatic  scattering  coefficient  of  random  rough  surfaces  is  examined  in  this  paper.  The 
iterative  PO  method  offers  decreased  memory  and  computation  time  restrictions  compared  to  the  standard 
numerical  methods  such  as  the  Method  of  Moments  (MoM).  Results  from  the  iterative  PO  method  are 
compared  to  the  standard  electric  field  integral  equation  (EFIE),  the  magnetic  field  integral  equation 
(MFIE)  as  well  as  the  existing  theoretical  solutions  for  rough  surfaces.  It  is  demonstrated  that  memory 
requirements  and  computation  time  is  significantly  decreased  while  providing  fairly  accurate  results  for 
surfaces  with  moderate  to  low  rms  slope. 

I  Introduction 

Development  of  numerically  efficient  Monte  Carlo-based  models  for  simulations  of  electromagnetic  scattering 
from  random  surfaces  has  attained  significant  prominence  over  the  past  decade  [4,  2,  1].  A  major  stumbling 
block  in  this  endeavor  has  been  the  large  memory  and  computation  time  requirement.  This  is  because  the  size 
of  the  scatter  is  large  compared  to  the  wavelength  and  the  Monte  Carlo  simulations  demand  computation  of 
the  scattering  problem  many  times.  Iterative  methods  offer  an  alternative  approach  when  exact  solutions  are 
not  available  and  have  been  used  in  different  electromagnetic  problems.  By  construct  evaluation  of  iterative 
solutions  are  rather  straight  forward  especially  when  the  perturbation  parameter  is  relatively  small.  Physical 
Optics  (PO)  approximation  is  known  to  provide  accurate  approximation  for  the  induced  surface  currents 
provided  that  the  local  radii  of  curvature  at  each  point  on  the  surface  of  scatterer  is  large  and  the  surface  is 
convex.  For  concave  surfaces  and  surfaces  with  many  adjacent  humps  multiple  scattering  drastically  alter  the 
standard  PO  current.  However  these  surface  current  variations  can  be  estimated  through  an  iterative  process. 
Unlike  Method  of  Moments  (MoM)  which  requires  matrices  on  the  order  of  A2  to  find  the  surface  current  of 
the  sample  surface  with  N  elements,  the  iterative  Physical  Optics  method  only  requires  memory  size  of  the 
order  N.  Thus,  substantial  memory  savings  are  realized.  Also,  since  no  solver  routine  is  necessary  in  order  to 
solve  for  the  surface  currents,  as  in  the  MoM,  substantial  time  savings  are  realized  as  well. 

We  investigated  the  use  of  the  iterative  Physical  Optics  method  upon  a  variety  of  surfaces  in  order  to  find  the 
approximate  region  of  validity  for  such  a  method.  The  results  obtained  using  the  iterative  PO  method  are  also 
compared  to  the  results  found  using  the  electric  field  integral  equation  (EFIE)  with  tapered  resistive  sheets 
at  the  ends  of  the  surface  samples  as  well  as  the  magnetic  field  integral  equation  (MFIE)  for  a  horizontally 
polarized  wave. 

II  Formulation 

Using  a  Monte  Carlo  simulation,  first  many  sample  surfaces  of  a  stochastic  random  process  representing  the 
desired  random  rough  surface  is  generated  using  a  standard  procedure  [4,  2].  As  mentioned  earlier  it  is  far 
more  convenient  to  solve  for  the  induced  polarization  current  using  an  iterative  method  as  opposed  to  brute 
force  numerical  methods.  Magnetic  field  integral  equation  (MFIE)  can  be  used  as  the  basis  for  the  iterative 
PO  solution  using  the  surface  curvature  as  the  perturbation  parameters.  The  MFIE  for  the  induced  surface 
current  (J3)  over  a  perfect  electric  conductor  is  given  by  [6] 


J,(r)  =  2  <»  x  H<)  +  i  /j,(r')(»  •  (r-  -  r»  (it  -  ^L.) 


(1) 


where  H1  is  the  incident  magnetic  field,  h  is  the  unit  normal  at  the  observation  point,  and  f  represents 
the  principle  value  integral.  It  is  obvious  from  (1)  that  for  a  flat  surface  (r'  -  r)  •  n  =  0  which  renders 
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Js(r)  —  2 (n  x  Hl)  which  is  the  standard  Physical  Optics  approximation  for  the  induced  surface  current.  For 
undulating  surfaces  with  large  radii  of  curvature  the  contribution  from  the  integral  has  a  secondary  effect  and 
therefore  (1)  may  be  solved  in  an  in  iterative  fashion.  To  examine  the  accuracy  of  the  iterative  PO  approach  the 
results  must  be  compared  with  an  exact  solution  obtained  from  a  numerical  method.  Numerical  solutions  for 
one-dimensional  rough  surfaces  are  tractable  and  thus  the  accuracy  of  the  iterative  PO  approach  is  examined 
for  one-dimensional  rough  surfaces.  The  MFIE  for  two-dimensional  problems  (one-dimensional  roughness)  can 
easily  be  obtained.  For  2-D  scattering  problems  the  transverse  ( Jt )  and  longitudinal  (Jz)  components  of  the 
induced  current  are  decoupled  and  the  integral  equations  for  a  TM  and  TE  incidence  are  respectively  given 
by: 


Mp)  =  2(n  x  H  i)-z+~  fMp')H[l\k0\p-J\)n'{*_  ^  dt!  (2) 

-MP)  =  2 +  y-  fjt{p')H?\h\p  -  dC>  (3) 

where  t  =  z  x  h  is  the  tangent  unit  vector.  Again  we  recognize  that  for  a  flat  surface  (h  =  z)  the  contribution 
from  the  integrals  in  (2)  and  (3)  are  zero  and  the  PO  currents  are  retrieved. 

Consider  a  perfectly  conducting  rough  surface  illuminated  by  an  E-polarized  (TM)  plane  wave.  Points  on  the 
surface  can  be  grouped  into  two  categories:  1)  lit  points  and  2)  shadowed  points.  As  a  first-order  approximation 
the  induced  current  over  the  lit  region  is  given  by  ji1^  =  2(h  x  H*)  and  over  the  shadow  regions  =  0.  This 
current  produces  a  scattered  magnetic  field  whose  tangent  on  the  surface  induces  a  secondary  PO  current.  It 
can  be  shown  that  the  expression  for  this  first  order  scattered  field  is  given  by 

*  *  H(i)  =  x  fj?\p')Hi\**\p-p'\)- dt  (4) 

and  hence  the  second  order  PO  current  can  be  obtained  from  J ^  =  2 (h  x  H^)  which  is  exactly  the  same  as 

the  integral  in  (2).  From  this  argument  it  is  evident  that  starting  with  Jzl\  the  integral  in  (2)  can  be  used  in 
a  recursive  manner  to  find  a  higher  order  solution.  There  is  a  subtlety  in  the  computation  of  the  second  order 
current  over  the  shadow  regions.  The  reason  for  shadow  is  that  the  scattered  field  is  out  of  phase  with  the 
incident  wave.  Basically  we  have  to  add  the  contribution  from  the  incident  field  to  the  first-order  scattered 
fields  over  the  shadow  regions.  Therefore  the  corrected  second  order  solution  is 

j(2)corrected  _  f  over  lit  region  ^ 

1  ,/i2'1  +  2  (n  x  H*)  over  shadow  region  v 

characterized,  higher  order  currents  can  be  obtained  from 

j(n)  =  jj(n- l)ip>)HW{kolp_p>l)lt_k^l  d(!  (6) 

The  convergence  is  examined  by  monitoring  the  normalized  difference  in  the  successive 


III  MFIE  and  the  Method  of  Moments 

Using  Galerkin’s  Method,  equation  (2)  is  re-written  as 

N 

j*  ^m2[h  x  H’(p)]  •  zdt  =  ^ <j>mJz(p)  —  ^  j^j>m4>nJz{p')H^\ko\p  —  d£d£'  .  (7) 


Once  jl2)corrected  is 


and  Jz  =  J2n 
(\ 
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which  can  be  solved  using  numerical  integration,  and  cast  into  a  matrix  equation.  In  equation  (7),  on  and  om 
are  the  basis  function  and  testing  function  respectively. 

This  is  also  compared  with  the  EFIE  which,  using  Galerkin's  method  is  given  by 

j[w (p)dc  =  R{p)3e{p)  +  ^^J^m4>^^p)HiY\k0\p-p'\)dCdC  (8) 

where  R  is  the  surface  resistivity  (for  a  perfect  conductor,  R  =  0).  For  the  EFIE,  tapered  resistive  sheets  are 
used  to  reduce  the  edge  effects.  An  optimum  tapered  function  is  reported  in  [2]. 

IV  Results 

To  verify  the  accuracy  of  the  iterative  PO  method,  a  perfectly  conducting  surface  with  two  Gaussian  humps,  as 
shown  in  Figure  1(a),  was  studied.  The  total  bistatic  scattering  pattern  for  normal  incidence  excitation  is  shown 
in  Figure  1(b)  and  excellent  agreement  is  observed  between  the  electric  field  integral  equation  (EFIE),  iifagne tic 
field  integral  equation  (MFIE)  and  the  iterative  Physical  Optics  (PO)  method.  After  verifying  that  the  results 
from  the  iterative  PO  method,  EFIE  and  MFIE  gave  similar  results,  numerical  Monte  Carlo  simulations 
of  rough  surface  scattering  is  performed  for  four  different  surfaces  denoted  by  Sx,  S2,  £3,  and  S4.  The 
roughness  parameters  ( ks ,  k£)  for  each  of  these  surfaces,  created  using  [4],  are  respectively  k$\  —  0.3,  k£\  =  3.0; 
ks2  =  0.5,  H2  =  6.13;  ks3  =  1.0,  k£$  =  3.0;  ks±  =  0.5,  &I4  =  2.0.  Surfaces  S\  and  S2  were  chosen  because 
they  fall  into  the  regions  of  validity  for  the  Small  Perturbation  Method  (SPM)  and  the  Physical  Optics  (PO) 
method,  respectively. 

For  surfaces  Si  and  S2,  whose  bistatic  scattering  patterns  (cr°)  are  shown  in  Figures  2(a)  and  2(b)  respectively, 
good  agreement  is  shown  between  the  EFIE,  MFIE,  iterative  PO  method  and  the  analytical  solutions.  Si 
sample  surfaces  are  30  A  in  length  (except  for  the  EFIE  where  1  A  resistive  tapered  sheets  are  added  on  each 
end)  with  Ax  =  0.2A  and  0,-  =  30.0°.  40  sample  surfaces  are  used  to  find  the  estimate  of  <r°.  Good  agreement 
is  shown  in  the  figure  until  about  -60°.  S2  sample  surfaces  are  46  A  in  length  with  Ax  =  0.1  A  and  0,  =  30.0°. 
Again,  good  agreement  is  shown  in  the  figure  until  about  -60°. 

For  surfaces  S3  and  S4,  good  agreement  is  also  shown  for  the  bistatic  scattering  pattern.  For  surface  S3, 
N  =  40,  length  =  18A,  0,-  =  30.0°,  rms  slope  =  0.471  and  Ax  =  0.1A,  and  the  bistatic  scattering  pattern  is 
shown  in  Figure  3(a).  For  surface  S4j  N  =  40, length  =  18A,  0,-  =  30.0°,  rms  slope  =  0.354  and  Ax  =  0.1A, 
and  the  bistatic  scattering  pattern  is  shown  in  Figure  3(b). 

V  Conclusion 

It  has  been  found  that  using  an  iterative  PO  method  for  characterizing  the  bistatic  scattering  pattern  for 
random  rough  surfaces  with  low  to  moderate  rms  slope  (<  0.6),  relatively  good  agreement  with  the  exact 
solution  using  both  the  EFIE  and  the  MFIE  is  demonstrated.  Computation  time  is  significantly  decreased 
as  well  as  memory  size  necessary  to  perform  the  numerical  solution  for  the  iterative  PO  method.  Since  for 
3-D  problems  the  kernel  of  the  integral  equation  decays  faster  with  distance  than  that  for  2-D  problems,  it  is 
expected  that  iterative  P.O.  to  perform  even  better  for  3-D  problems. 
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Figure  1:  (a)  Surface  profile  with  two  Gaussian  humps,  (b)  Total  bistatic  scattering 


(a)  (b) 

Figure  2:  (a)  Comparison  of  bistatic  scattering  from  a  random  surface  with  ks  =  0.3  kt  —  3.0,  N  =  40, 
length  =  30A,  resistive  tapered  ends  length  for  EFIE  =  1A,  Ax  =  0.2A  and  0,-  =  30°.  Comparison  shows 
difference  between  EFIE,  MFIE,  iterative  PO  method,  and  the  analytical  solution,  (b)  Comparison  of  bistatic 
scattering  from  a  random  surface  ks  =  0.5  kt  —  6.13,  N  —  40,  length  =  46A,  resistive  tapered  ends  length  for 
EFIE  =  1A,  Ax  =  0.1  A  and  0*  =  30°.  Comparison  shows  difference  between  EFIE,  iterative  PO  method,  and 
the  Kirchoff  Approximation. 


(a)  (b) 

Figure  3:  (a)  Comparison  of  bistatic  scattering  from  a  random  surface  with  ks  =  1.0  kt  =  3.0,  N  =  40, 
length  =  18A,  resistive  tapered  ends  length  for  EFIE  =  1A,  Ax  —  0.1  A  and  0,-  =  30.0°.  Comparison  shows 
difference  between  EFIE,  MFIE,  and  the  iterative  PO  method,  (b)  Comparison  of  bistatic  scattering  from  a 
random  surface  with  ks  =  0.5  kt  =  2.0,  N  =  40,  length  =  18A,  resistive  tapered  ends  length  for  EFIE  =  1A, 
Ax  =  0.1A  and  0;  =  30.0°.  Comparison  shows  difference  between  EFIE,  MFIE,  and  the  iterative  PO  method. 
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Abstract 

The  application  of  a  fast  far-field  iterative  Physical  Optics  (FIPO)  method  in  conjunction  with  a  Monte 
Carlo  simulation  for  characterizing  the  bistatic  scattering  coefficient  of  random  rough  surfaces  is  examined 
in  this  paper.  The  FIPO  method  offers  decreased  memory  and  computation  time  restrictions  compared  to 
the  standard  numerical  methods  such  as  the  Method  of  Moments  (MoM),  and  decreased  computation  time 
compared  to  an  exact  iterative  PO  method.  Results  from  the  FIPO  method  are  compared  to  the  standard 
electric  field  integral  equation  (EFIE),  the  magnetic  field  integral  equation  (MFIE),  a  complete  iterative 
PO  (IPO),  as  well  as  the  existing  theoretical  solutions  for  rough  surfaces.  It  is  demonstrated  that  memory 
requirements  and  computation  time  is  significantly  decreased  while  providing  fairly  accurate  results  for 
surfaces  with  moderate  to  low  rms  slope. 


I  Introduction 

Numerically  efficient  Monte  Carlo-based  models  for  simulations  of  electromagnetic  scattering  from  random 
surfaces  has  attained  significant  prominence  over  the  past  decade  [4,  2,  1].  A  major  stumbling  block  in  this 
endeavor  has  been  the  large  memory  and  computation  time  requirement,  but  as  previously  reported,  the 
computation  time  and  memory  size  has  been  reduced  using  a  standard  iterative  PO  method  which  produces 
fairly  accurate  results  [9].  For  a  surface  consisting  of  N  elements,  the  iterative  Physical  Optics  method  only 
requires  memory  size  of  the  order  N.  Thus,  substantial  time  and  memory  savings  are  realized  compared  to 
the  Method  of  Moments  (MoM). 

Previously,  we  investigated  the  use  of  the  iterative  Physical  Optics  method  upon  a  variety  of  surfaces  in 
order  to  find  the  approximate  region  of  validity  for  such  a  method  [9].  The  difference  between  a  complete 
iterative  PO  approach  and  the  new  fast  far-field  IPO  (FIPO)  approach  is  that  the  original  IPO  was  order 
N 2  in  computation  time,  while  the  new  FIPO  approach  is  only  order  N  computation  time.  Both  routines 
are  order  N  in  memory  requirements.  The  results  obtained  using  the  FIPO  method  are  also  compared  to  the 
results  found  using  the  electric  field  integral  equation  (EFIE)  with  tapered  resistive  sheets  at  the  ends  of  the 
surface  samples,  the  magnetic  field  integral  equation  (MFIE),  and  the  complete  IPO  method  for  a  horizontally 
polarized  wave. 


II  Formulation 


For  a  Monte  Carlo  simulation,  many  sample  surfaces  representing  the  desired  random  rough  surface  character¬ 
istics  are  generated  using  a  standard  surface  generation  routine  [4,  2].  As  previously  reported,  the  Magnetic 
field  integral  equation  (MFIE)  is  used  as  the  basis  for  the  iterative  PO  method  [9]  and  the  currents  over  a 
perfect  electric  conductor  is  given  by  [6] 


Js(r)  =  2  (»  x  H*)  +  T  fj,(v')(h  ■  (r'  -  r))  (ik  -  — 


gifcflr— r'|) 

|r-  r'l2 


ds' 


(1) 


where  H‘  is  the  incident  magnetic  field,  n  is  the  unit  normal  at  the  observation  point,  and  /  represents  the 
principle  value  integral.  For  the  MFIE  for  two-dimensional  scattering  problems  (one-dimensional  roughness), 
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the  transverse  ( Jt )  and  longitudinal  (Jz)  components  of  the  induced  current  are  decoupled  and  the  integral 
equations  for  a  TM  and  TE  incidence  are  respectively  given  by: 

Up)  =  2(n  X  H1)  •  i  +  ^  -  p'\)n  ^P_~p^]  dC  (2) 

-Jt(p)  =  2  Hi  +  '-j-  jjt(p')H[l)(k0\p  -  P'lf' dC  (3) 

where  t !  =  £  x  h  is  the  tangent  unit  vector. 

Consider  a  perfectly  conducting  rough  surface  illuminated  by  an  E- polarized  (TM)  plane  wave.  Points  on  the 
surface  are  grouped  into  two  categories:  1)  lit  points  and  2)  shadowed  points.  As  a  first-order  approximation 
the  induced  current  over  the  lit  region  is  given  by  =  2(n  x  Hf)  and  over  the  shadow  regions  =  0.  This 
current  produces  a  scattered  magnetic  field  whose  tangent  on  the  surface  induces  a  secondary  PO  current. 
The  expression  for  this  first  order  scattered  field  is  given  by  - . 

»  x  Hfu  =  ^  fjV{p')H?\k0\p  -  p'\) n dt  (4) 

and  hence  the  second  order  PO  current  can  be  obtained  from  J j2^  =  2(h  x  H^)  which  is  exactly  the  same  as 
the  integral  in  (2). 

For  the  fast  far-field  iterative  PO  (FIPO),  instead  of  a  complete  solution  as  given  in  (2)  or  (3),  the  surface  is 
divided  into  a  near-field  region  and  far-field  regions.  Then,  (2)  and  (3)  can  be  written  as 


Ji ">(»)  =  [  j(--)(p')gW(t0|),-(,'|)^7i<’-r  p'>  de 

Jnearfteldm  \P  P\  ^  Jfarfieldmt  \P  P\ 


(5) 


=  \t,-P'\)^lELjzAde+  £  f 

J  nearfieldm  \P-P'\  m=1  .  2  Jfarfieldm,  I P~ 


~p) 

P'1 

(6) 


Then,  the  FIPO  can  be  found  by  summing  the  contributions  from  the  near-field  plus  the  far-field  contributions. 

The  far-field  contributions,  or  last  terms,  from  (5)  and  (6)  may  be  approximated  respectively  using  the  large 

argument  expansion  of  Hankel  functions  and  noting  that  in  the  far  field  ,p~pti  &  u  where  u  =  x  if  o'  >  p~  or 
„  I p-p  I  * 

u  =  —x  if  <  px 


-  P'mid |)  /  Jrt(B"1)(p0e-,*°^»“-'')- V  •  *)<&  (7) 

Z  J  far  field 

^H^ikolp  -  p'mid\)(h  •  u)  f  j(n-'){p')e-iko^<-e‘y*dt'  (8) 

Z  J  jar  field 

where  Pmid  is  the  midpoint  of  the  far  field  section.  For  the  far  field,  since  the  fields  fall  off  quickly  due  to 
the  kernel,  large  sections  of  surface  are  integrated  over  to  give  us  the  far-field  contribution  to  the  integral.  As 
shown  in  (7),  for  any  near-field  section,  the  far-field  section  only  has  to  be  calculated  once.  Then  instead  of 
an  N 2  calculation  the  order  of  calculation  is  NM ,  that  is,  a  computational  savings  factor  of  (-“)  is  achieved, 
where  M  is  the  total  number  of  far-field  sections.  As  an  example,  for  a  40A  surface  consisting  of  400  elements 
divided  every  5A  which  leads  to  m  =  8  distinct  sections,  FIPO  is  approximately  50  times  faster  than  the 
traditional  IPO.  In  the  actual  implementation,  the  adjacent  sections  should  also  be  excluded  from  the  far-field 


■dt 
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approximation,  that  is,  mf  ^  {w  —  1.777.77?  +  1}.  Then,  for  the  example  above,  for  each  section  ???,  there  are 
6  far- field  sections  if  m  is  equal  to  1  or  8,  and  5  far-field  sections  if  in  is  in  the  range  {2.. .7}.  For  any  section, 
m,  there  are  100  or  150  near-field  elements  which  have  to  be  iterated  over,  and  5  or  6  far-field  sections,  which 
are  calculated  once. 

The  above  argument  allows  one  to  start  with  and  recursively  find  a  higher  order  solution.  As  reported 
previously  there  exists  a  subtlety  in  the  computation  of  the  second  order  current  over  the  shadow  regions  [9]. 
The  corrected  second  order  solution  is 


j(2)corrected  _ 


/  o\ 

Jz  over  lit  region 

A ^  +  2  (7?  x  H1)  over  shadow  region 


(9) 


From  j(2)corrected }  higher  order  currents  can  be  obtained. 

The  convergence  still  is  examined  by  monitoring  the  normalized  difference  in  the  successive  solutions 
as  in  [9]. 


( 


£n=, 


III  Results 

Previously,  it  was  determined  that  the  iterative  PO  method  produced  very  favorable  results  for  surfaces  with 
small  RMS  slope.  The  FIPO  method  was  then  compared  to  the  IPO  method  for  a  number  of  surfaces  and 
found  to  also  produce  very  favorable  result.  These  results  are  compared  to  the  complete  IPO,  the  Electric 
Field  Integral  Equation  (EFIE)  as  well  as  the  Magnetic  Field  Integral  Equation  (MFIE). 

First,  the  FIPO  method  was  compared  to  the  EFIE,  MFIE,  and  IPO  method  for  a  single  Gaussian  hump, 
as  shown  in  Figure  la.  The  bistatic  scattering  pattern  for  a  single  Gaussian  hump  is  shown  in  Figure  lb 
and  excellent  agreement  between  the  EFIE,  MFIE,  IPO,  and  FIPO  method  are  shown.  After  verifying  these 
results,  two  surfaces,  Si  and  S2  were  compared  for  the  IPO  method  as  well  as  the  FIPO  method.  Surface  S\ 
is  a  surface  that  was  measured  in  [9]  and  produced  excellent  agreement  with  the  EFIE  and  MFIE  method.  Si 
has  surface  characteristics  ks  =  0.5  and  kt  =  2.0.  In  addition,  a  Monte  Carlo  simulation  of  surfaces  that  have 
characteristics  ks  =  0.65,  kt  =  2.0  and  where  each  surface  is  800A  in  length  (or  8000  elements)  is  compared. 

For  surface  SI,  the  total  bistatic  scattering  pattern  (<7o)  is  shown  in  Figure  2a.  S\  sample  surfaces  are  18A 
in  length,  with  Ax  =  0.1  A,  the  number  of  surfaces  averaged  over,  N ,  is  40  and  0*=3O.O°.  For  surface  S2 ,  the 
total  bistatic  scattering  pattern  is  shown  in  Figure  2b,  and  the  IPO  and  FIPO  methods  are  compared  to  result 
obtained  by  surfaces  with  the  same  characteristics  ( ks ,  kt).  It  is  obvious  the  both  the  IPO  and  FIPO  methods 
work  very  well  compared  to  the  EFIE  and  MFIE  methods.  For  surface  S2}  ks  =  0.65,  kt  =  2.0,  N  =  30, length 
=  800A  for  FIPO  and  IPO,  while  length  =  20A  for  EFIE,  and  =  30.0°. 


IV  Conclusion 


It  has  been  found  that  when  using  a  fast  far-field  iterative  PO  method  for  characterizing  the  bistatic  scattering 
pattern  for  random  rough  surfaces  with  low  to  moderate  rms  slope  (<  0.6),  relatively  good  agreement  with 
the  complete  iterative  PO  method,  as  well  as  with  the  exact  solutions  using  both  the  EFIE  and  the  MFIE  is 
demonstrated.  Computation  time  is  significantly  decreased  as  well  as  memory  size  necessary  to  perform  the 
numerical  solution  for  the  iterative  PO  method.  Since  for  3-D  problems  the  kernel  of  the  integral  equation 
decays  faster  with  distance  than  that  for  2-D  problems,  it  is  expected  that  iterative  PO  to  perform  even  better 
for  3-D  problems. 
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Figure  2:  (a)  Comparison  of  bistatic  scattering  from  a  random  surface  with  ks  =  1.0  kt  =  3.0,  N  =  40,  length 
=  18A,  resistive  tapered  ends  length  for  EFIE  =  1A,  Ax  =  0.1A  and  0,-  =  30.0°.  Comparison  shows  difference 
between  IPO,  FIPO,  and  EFIE  method,  (b)  Comparison  of  bistatic  scattering  from  a  random  surface  using 
FIPO  and  IPO  ks  =  0.65  ,kl  =  2.0,  N  =  20,  length  =  800A,  Ax  =  0.1A  and  0,-  =  30.0°  vs.  the  EFIE  for 
the  same  physical  characteristics,  but  IV  =  40,  length  =  20 A  with  tapered  resistive  sheets  of  1A.  Comparison 
shows  difference  between  IPO,  FIPO,  and  EFIE  method. 
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